+Initialize thicknesses in height units

  Pass arguments in height units rather than thickness units to most of the
routines that initialize thickness or temperatures and salinities.  These
routines are already undoing this scaling and working in height units, and it is
not possible to convert thicknesses to thickness units in non-Boussinesq mode
until the temperatures and salinities are also known.  The routines whose
argument units are altered include:

- initialize_thickness_uniform
- initialize_thickness_list
- DOME_initialize_thickness
- ISOMIP_initialize_thickness
- benchmark_initialize_thickness
- Neverworld_initialize_thickness
- circle_obcs_initialize_thickness
- lock_exchange_initialize_thickness
- external_gwave_initialize_thickness
- DOME2d_initialize_thickness
- adjustment_initialize_thickness
- sloshing_initialize_thickness
- seamount_initialize_thickness
- dumbbell_initialize_thickness
- soliton_initialize_thickness
- Phillips_initialize_thickness
- Rossby_front_initialize_thickness
- user_initialize_thickness
- DOME2d_initialize_temperature_salinity
- ISOMIP_initialize_temperature_salinity
- adjustment_initialize_temperature_salinity
- baroclinic_zone_init_temperature_salinity
- sloshing_initialize_temperature_salinity
- seamount_initialize_temperature_salinity
- dumbbell_initialize_temperature_salinity
- Rossby_front_initialize_temperature_salinity
- SCM_CVMix_tests_TS_init
- dense_water_initialize_TS
- adjustEtaToFitBathymetry

Similar changes were made internally to MOM_temp_salt_initialize_from_Z to defer
the transition to working in thickness units, although the appropriate call to
convert_thickness does still occur within MOM_temp_salt_initialize_from_Z and
the units of its arguments are not changed.

  The routine convert thickness was modified to work with a new input depth
space input thickness argument and return a thickness in thickness units, and it
is now being called after all of the routines to initialize thicknesses and
temperatures and salinities, except in the few cases where the thickness are
being specified directly in mass-based thickness units, as might happen when
they are read from an input file.

  The new option "mass_file" is now a recognized option for the THICKNESS_CONFIG
runtime parameter, and this information is passed in the new mass_file argument
to initialize_thickness_from_file.  The description of the runtime parameter
THICKNESS_IC_RESCALE was updated to reflect this change.

  The unused thickness (h) argument to soliton_initialize_velocity was
eliminated.

  The unused thickness (h) argument to determine_temperature was eliminated, as
was the unused optional h_massless argument to the same function.

  This commit also rearranges the calls to do adjustments to the thicknesses to
account for the presence of an ice shelf or to iteratively apply the ALE
remapping to occur before the velocities are initialized, so that there is a
clearer separation of the phases of the initialization.

  Also added optional height_units argument to ALE_initThicknessToCoord to
specify that the coordinate are to be returned in height_units.  If it is
omitted or false, the previous thickness units are returned, but when called
from MOM_initialize_state the new argument is being used.

  The runtime parameter CONVERT_THICKNESS_UNITS is no longer meaningful, so it
has been obsoleted.

  All answers are bitwise identical, but there are multiple changes to the
arguments to publicly visible subroutines or their units, and there are changes
to the contents of the MOM_parameter_doc files.

src/ALE/MOM_ALE.F90

@@ -1463,17 +1463,23 @@ subroutine ALE_writeCoordinateFile( CS, GV, directory )
 end subroutine ALE_writeCoordinateFile
 
 !> Set h to coordinate values for fixed coordinate systems
-subroutine ALE_initThicknessToCoord( CS, G, GV, h )
+subroutine ALE_initThicknessToCoord( CS, G, GV, h, height_units )
   type(ALE_CS), intent(inout)                            :: CS  !< module control structure
   type(ocean_grid_type), intent(in)                      :: G   !< module grid structure
   type(verticalGrid_type), intent(in)                    :: GV  !< Ocean vertical grid structure
-  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h   !< layer thickness [H ~> m or kg m-2]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h   !< layer thickness in thickness units
+                                                                !! [H ~> m or kg m-2] or height units [Z ~> m]
+  logical,                          optional, intent(in) :: height_units !< If present and true, the
+                                                                !! thicknesses are in height units
 
   ! Local variables
+  real :: scale ! A scaling value for the thicknesses [nondim] or [H Z-1 ~> nondim or kg m-3]
   integer :: i, j
 
+  scale = GV%Z_to_H
+  if (present(height_units)) then ; if (height_units) scale = 1.0 ; endif
   do j = G%jsd,G%jed ; do i = G%isd,G%ied
-    h(i,j,:) = GV%Z_to_H * getStaticThickness( CS%regridCS, 0., G%bathyT(i,j)+G%Z_ref )
+    h(i,j,:) = scale * getStaticThickness( CS%regridCS, 0., G%bathyT(i,j)+G%Z_ref )
   enddo ; enddo
 
 end subroutine ALE_initThicknessToCoord

src/diagnostics/MOM_obsolete_params.F90

@@ -56,6 +56,7 @@ subroutine find_obsolete_params(param_file)
          hint="Instead use OBC_SEGMENT_xxx_VELOCITY_NUDGING_TIMESCALES.")
   enddo
 
+  call obsolete_logical(param_file, "CONVERT_THICKNESS_UNITS", .true.)
   call obsolete_logical(param_file, "MASK_MASSLESS_TRACERS", .false.)
 
   call obsolete_logical(param_file, "SALT_REJECT_BELOW_ML", .false.)

src/tracer/MOM_tracer_Z_init.F90

@@ -556,29 +556,24 @@ end function find_limited_slope
 
 !> This subroutine determines the potential temperature and salinity that
 !! is consistent with the target density using provided initial guess
-subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, h, k_start, G, GV, US, &
-                                 PF, just_read, h_massless)
+subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, k_start, G, GV, US, PF, &
+                                 just_read)
   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(inout) :: temp !< potential temperature [C ~> degC]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                                  intent(inout) :: salt !< salinity [S ~> ppt]
   real, dimension(SZK_(GV)),     intent(in)    :: R_tgt !< desired potential density [R ~> kg m-3].
-  type(EOS_type),                intent(in)    :: EOS !< seawater equation of state control structure
+  type(EOS_type),                intent(in)    :: EOS   !< seawater equation of state control structure
   real,                          intent(in)    :: p_ref !< reference pressure [R L2 T-2 ~> Pa].
   integer,                       intent(in)    :: niter !< maximum number of iterations
   integer,                       intent(in)    :: k_start !< starting index (i.e. below the buffer layer)
-  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
-                                 intent(in)    :: h   !< layer thickness, used only to avoid working on
-                                                      !! massless layers [H ~> m or kg m-2]
   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.
   logical,                       intent(in)    :: just_read !< If true, this call will only read
                                                       !! parameters without changing T or S.
-  real,                optional, intent(in)    :: h_massless !< A threshold below which a layer is
-                                                      !! determined to be massless [H ~> m or kg m-2]
 
   ! Local variables (All of which need documentation!)
   real, dimension(SZI_(G),SZK_(GV)) :: &
@@ -587,7 +582,6 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, h, k_star
     dT, &  ! An estimated change in temperature before bounding [C ~> degC]
     dS, &  ! An estimated change in salinity before bounding [S ~> ppt]
     rho, & ! Layer densities with the current estimate of temperature and salinity [R ~> kg m-3]
-    hin, & ! A 2D copy of the layer thicknesses [H ~> m or kg m-2]
     drho_dT, & ! Partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1]
     drho_dS    ! Partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]
   real, dimension(SZI_(G)) :: press ! Reference pressures [R L2 T-2 ~> Pa]
@@ -675,7 +669,6 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, h, k_star
     dS(:,:) = 0. ! Needs to be zero everywhere since there is a maxval(abs(dS)) later...
     T(:,:) = temp(:,j,:)
     S(:,:) = salt(:,j,:)
-    hin(:,:) = h(:,j,:)
     dT(:,:) = 0.0
     adjust_salt = .true.
     iter_loop: do itt = 1,niter
@@ -685,7 +678,7 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, h, k_star
                                       EOS, EOSdom )
       enddo
       do k=k_start,nz ; do i=is,ie
-!       if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. hin(i,k)>h_massless .and. abs(T(i,k)-land_fill) < epsln) then
+!       if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. abs(T(i,k)-land_fill) < epsln) then
         if (abs(rho(i,k)-R_tgt(k))>tol_rho) then
           if (.not.fit_together) then
             dT(i,k) = max(min((R_tgt(k)-rho(i,k)) / drho_dT(i,k), max_t_adj), -max_t_adj)
@@ -713,7 +706,7 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, h, k_star
                                       EOS, EOSdom )
       enddo
       do k=k_start,nz ; do i=is,ie
-!       if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. hin(i,k)>h_massless .and. abs(T(i,k)-land_fill) < epsln ) then
+!       if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. abs(T(i,k)-land_fill) < epsln ) then
         if (abs(rho(i,k)-R_tgt(k)) > tol_rho) then
           dS(i,k) = max(min((R_tgt(k)-rho(i,k)) / drho_dS(i,k), max_s_adj), -max_s_adj)
           S(i,k) = max(min(S(i,k)+dS(i,k), S_max), S_min)

src/user/DOME2d_initialization.F90

@@ -98,7 +98,7 @@ subroutine DOME2d_initialize_thickness ( h, depth_tot, G, GV, US, param_file, ju
   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].
+                           intent(out) :: h           !< The thickness that is being initialized [Z ~> m]
   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
@@ -158,16 +158,16 @@ subroutine DOME2d_initialize_thickness ( h, depth_tot, G, GV, US, param_file, ju
           eta1D(k) = e0(k)
           if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then
             eta1D(k) = eta1D(k+1) + GV%Angstrom_Z
-            h(i,j,k) = GV%Angstrom_H
+            h(i,j,k) = GV%Angstrom_Z
           else
-            h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1))
+            h(i,j,k) = eta1D(k) - eta1D(k+1)
           endif
         enddo
 
         x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon
         if ( x <= dome2d_width_bay ) then
-          h(i,j,1:nz-1) = GV%Angstrom_H
-          h(i,j,nz) = GV%Z_to_H * dome2d_depth_bay * G%max_depth - (nz-1) * GV%Angstrom_H
+          h(i,j,1:nz-1) = GV%Angstrom_Z
+          h(i,j,nz) = dome2d_depth_bay * G%max_depth - (nz-1) * GV%Angstrom_Z
         endif
 
       enddo ; enddo
@@ -180,16 +180,16 @@ subroutine DOME2d_initialize_thickness ( h, depth_tot, G, GV, US, param_file, ju
  !         eta1D(k) = e0(k)
  !         if (eta1D(k) < (eta1D(k+1) + min_thickness)) then
  !           eta1D(k) = eta1D(k+1) + min_thickness
- !           h(i,j,k) = GV%Z_to_H * min_thickness
+ !           h(i,j,k) = min_thickness
  !         else
- !           h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1))
+ !           h(i,j,k) = eta1D(k) - eta1D(k+1)
  !         endif
  !       enddo
  !
  !       x = G%geoLonT(i,j) / G%len_lon
  !       if ( x <= dome2d_width_bay ) then
- !         h(i,j,1:nz-1) = GV%Z_to_H * min_thickness
- !         h(i,j,nz) = GV%Z_to_H * (dome2d_depth_bay * G%max_depth - (nz-1) * min_thickness)
+ !         h(i,j,1:nz-1) = min_thickness
+ !         h(i,j,nz) = dome2d_depth_bay * G%max_depth - (nz-1) * min_thickness
  !       endif
  !
  !    enddo ; enddo
@@ -202,16 +202,16 @@ subroutine DOME2d_initialize_thickness ( h, depth_tot, G, GV, US, param_file, ju
           eta1D(k) = e0(k)
           if (eta1D(k) < (eta1D(k+1) + min_thickness)) then
             eta1D(k) = eta1D(k+1) + min_thickness
-            h(i,j,k) = GV%Z_to_H * min_thickness
+            h(i,j,k) = min_thickness
           else
-            h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1))
+            h(i,j,k) = eta1D(k) - eta1D(k+1)
           endif
         enddo
       enddo ; enddo
 
     case ( REGRIDDING_SIGMA )
       do j=js,je ; do i=is,ie
-        h(i,j,:) = GV%Z_to_H*depth_tot(i,j) / nz
+        h(i,j,:) = depth_tot(i,j) / nz
       enddo ; enddo
 
     case default
@@ -225,11 +225,11 @@ end subroutine DOME2d_initialize_thickness
 
 !> Initialize temperature and salinity in the 2d DOME configuration
 subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, US, param_file, just_read)
-  type(ocean_grid_type),                     intent(in)  :: G !< Ocean grid structure
+  type(ocean_grid_type),                     intent(in)  :: G  !< Ocean 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 !< Potential temperature [C ~> degC]
-  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< Salinity [S ~> 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),SZK_(GV)), intent(out) :: T  !< Potential temperature [C ~> degC]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S  !< Salinity [S ~> ppt]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)  :: h  !< Layer thickness [Z ~> m]
   type(unit_scale_type),                     intent(in)  :: US !< A dimensional unit scaling type
   type(param_file_type),                     intent(in)  :: param_file !< Parameter file structure
   logical,                                   intent(in)  :: just_read !< If true, this call will
@@ -287,7 +287,7 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi
       do j=js,je ; do i=is,ie
         xi0 = 0.0
         do k = 1,nz
-          xi1 = xi0 + (GV%H_to_Z * h(i,j,k)) / G%max_depth
+          xi1 = xi0 + h(i,j,k) / G%max_depth
           S(i,j,k) = S_surf + 0.5 * S_range * (xi0 + xi1)
           xi0 = xi1
         enddo
@@ -298,7 +298,7 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi
       do j=js,je ; do i=is,ie
         xi0 = 0.0
         do k = 1,nz
-          xi1 = xi0 + (GV%H_to_Z * h(i,j,k)) / G%max_depth
+          xi1 = xi0 + h(i,j,k) / G%max_depth
           S(i,j,k) = S_surf + 0.5 * S_range * (xi0 + xi1)
           xi0 = xi1
         enddo

src/user/DOME_initialization.F90

@@ -105,7 +105,7 @@ subroutine DOME_initialize_thickness(h, depth_tot, G, GV, param_file, just_read)
   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].
+                           intent(out) :: h           !< The thickness that is being initialized [Z ~> m]
   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
@@ -141,9 +141,9 @@ subroutine DOME_initialize_thickness(h, depth_tot, G, GV, param_file, just_read)
       eta1D(K) = e0(K)
       if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then
         eta1D(K) = eta1D(K+1) + GV%Angstrom_Z
-        h(i,j,k) = GV%Angstrom_H
+        h(i,j,k) = GV%Angstrom_Z
       else
-        h(i,j,k) = GV%Z_to_H * (eta1D(K) - eta1D(K+1))
+        h(i,j,k) = eta1D(K) - eta1D(K+1)
       endif
     enddo
   enddo ; enddo

src/user/ISOMIP_initialization.F90

@@ -143,7 +143,7 @@ subroutine ISOMIP_initialize_thickness ( h, depth_tot, G, GV, US, param_file, tv
   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].
+                           intent(out) :: h           !< The thickness that is being initialized [Z ~> m]
   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
@@ -170,7 +170,7 @@ subroutine ISOMIP_initialize_thickness ( h, depth_tot, G, GV, US, param_file, tv
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
 
   if (.not.just_read) &
-    call MOM_mesg("MOM_initialization.F90, initialize_thickness_uniform: setting thickness")
+    call MOM_mesg("ISOMIP_initialization.F90, ISOMIP_initialize_thickness: setting thickness")
 
   call get_param(param_file, mdl,"MIN_THICKNESS", min_thickness, &
                  'Minimum layer thickness', units='m', default=1.e-3, do_not_log=just_read, scale=US%m_to_Z)
@@ -225,9 +225,9 @@ subroutine ISOMIP_initialize_thickness ( h, depth_tot, G, GV, US, param_file, tv
         eta1D(k) = e0(k)
         if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then
           eta1D(k) = eta1D(k+1) + GV%Angstrom_Z
-          h(i,j,k) = GV%Angstrom_H
+          h(i,j,k) = GV%Angstrom_Z
         else
-          h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1))
+          h(i,j,k) = eta1D(k) - eta1D(k+1)
         endif
       enddo
     enddo ; enddo
@@ -240,17 +240,17 @@ subroutine ISOMIP_initialize_thickness ( h, depth_tot, G, GV, US, param_file, tv
         eta1D(k) =  -G%max_depth * real(k-1) / real(nz)
         if (eta1D(k) < (eta1D(k+1) + min_thickness)) then
           eta1D(k) = eta1D(k+1) + min_thickness
-          h(i,j,k) = GV%Z_to_H * min_thickness
+          h(i,j,k) = min_thickness
         else
-          h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1))
+          h(i,j,k) = eta1D(k) - eta1D(k+1)
         endif
       enddo
     enddo ; enddo
 
   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 * depth_tot(i,j) / real(nz)
+      h(i,j,:) = depth_tot(i,j) / real(nz)
     enddo ; enddo
 
   case default
@@ -269,7 +269,7 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, depth_tot, G, GV, U
   type(unit_scale_type),                     intent(in)  :: US !< A dimensional unit scaling type
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T  !< Potential temperature [C ~> degC]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S  !< Salinity [S ~> 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),SZK_(GV)), intent(in)  :: h  !< Layer thickness [Z ~> m]
   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
@@ -334,10 +334,10 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, depth_tot, G, GV, U
       do j=js,je ; do i=is,ie
         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
+          xi0 = xi0 + 0.5 * h(i,j,k)     ! Depth in middle of layer
           S(i,j,k) = S_sur + dS_dz * xi0
           T(i,j,k) = T_sur + dT_dz * xi0
-          xi0 = xi0 + 0.5 * h(i,j,k) * GV%H_to_Z ! Depth at top of layer
+          xi0 = xi0 + 0.5 * h(i,j,k)     ! Depth at top of layer
         enddo
       enddo ; enddo
 
@@ -372,10 +372,10 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, depth_tot, G, GV, U
         xi0 = 0.0
         do k = 1,nz
           !T0(k) = T_Ref; S0(k) = S_Ref
-          xi1 = xi0 + 0.5 * h(i,j,k) * GV%H_to_Z
+          xi1 = xi0 + 0.5 * h(i,j,k)
           S0(k) = S_sur - dS_dz * xi1
           T0(k) = T_sur - dT_dz * xi1
-          xi0 = xi0 + h(i,j,k) * GV%H_to_Z
+          xi0 = xi0 + h(i,j,k)
           ! write(mesg,*) 'S,T,xi0,xi1,k',S0(k),T0(k),xi0,xi1,k
           ! call MOM_mesg(mesg,5)
         enddo
@@ -430,7 +430,7 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, depth_tot, G, GV, U
   !i=G%iec; j=G%jec
   !do k = 1,nz
   !  call calculate_density(T(i,j,k), S(i,j,k),0.0,rho_tmp,eqn_of_state, scale=US%kg_m3_to_R)
-  !  write(mesg,*) 'k,h,T,S,rho,Rlay',k,h(i,j,k),US%C_to_degC*T(i,j,k),US%S_to_ppt*S(i,j,k),rho_tmp,GV%Rlay(k)
+  !  write(mesg,*) 'k,h,T,S,rho,Rlay',k,US%Z_to_m*h(i,j,k),US%C_to_degC*T(i,j,k),US%S_to_ppt*S(i,j,k),rho_tmp,GV%Rlay(k)
   !  call MOM_mesg(mesg,5)
   !enddo