+Add grid_unit_to_L to the ocean_grid_type

  Add the new element grid_unit_to_L to the ocean_grid_type and the
dyn_horgrid_type, which can be used to convert the units of the geoLat and
geoLon variables to rescaled horizontal distance units ([L ~> m]) when they are
Cartesian coordinates.  When Cartesian coordinates are not in use,
G%grid_unit_to_L is set to 0.

  This new element of the grid type is used to test for inconsistent grids or to
eliminate rescaling variables in set_rotation_beta_plane(),
initialize_velocity_circular(), DOME_initialize_topography(),
DOME_initialize_sponges(), DOME_set_OBC_data(), ISOMIP_initialize_topography(),
idealized_hurricane_wind_forcing(), Kelvin_set_OBC_data(),
Rossby_front_initialize_velocity(), soliton_initialize_thickness(), and
soliton_initialize_velocity().   These are the instances where this new
variable could be used and bitwise identical answers are recovered.  There are
a few other places where they should be used, but where answers would change,
and these will be deferred to a subsequent commit.

  All answers are bitwise identical, but there are new elements in two
transparent and widely used types.

src/core/MOM_grid.F90

@@ -190,6 +190,10 @@ module MOM_grid
 
   ! These parameters are run-time parameters that are used during some
   ! initialization routines (but not all)
+  real :: grid_unit_to_L !< A factor that converts a the geoLat and geoLon variables and related
+                        !! variables like len_lat and len_lon into rescaled horizontal distance
+                        !! units on a Cartesian grid, in [L km ~> 1000] or [L m-1 ~> 1] or
+                        !! is 0 for a non-Cartesian grid.
   real :: south_lat     !< The latitude (or y-coordinate) of the first v-line [degrees_N] or [km] or [m]
   real :: west_lon      !< The longitude (or x-coordinate) of the first u-line [degrees_E] or [km] or [m]
   real :: len_lat       !< The latitudinal (or y-coord) extent of physical domain [degrees_N] or [km] or [m]

src/core/MOM_transcribe_grid.F90

@@ -135,6 +135,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US)
   ! Copy various scalar variables and strings.
   oG%x_axis_units = dG%x_axis_units ; oG%y_axis_units = dG%y_axis_units
   oG%x_ax_unit_short = dG%x_ax_unit_short ; oG%y_ax_unit_short = dG%y_ax_unit_short
+  oG%grid_unit_to_L = dG%grid_unit_to_L
   oG%areaT_global = dG%areaT_global ; oG%IareaT_global = dG%IareaT_global
   oG%south_lat = dG%south_lat ; oG%west_lon  = dG%west_lon
   oG%len_lat = dG%len_lat ; oG%len_lon = dG%len_lon
@@ -296,6 +297,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US)
   ! Copy various scalar variables and strings.
   dG%x_axis_units = oG%x_axis_units ; dG%y_axis_units = oG%y_axis_units
   dG%x_ax_unit_short = oG%x_ax_unit_short ; dG%y_ax_unit_short = oG%y_ax_unit_short
+  dG%grid_unit_to_L = oG%grid_unit_to_L
   dG%areaT_global = oG%areaT_global ; dG%IareaT_global = oG%IareaT_global
   dG%south_lat = oG%south_lat ; dG%west_lon  = oG%west_lon
   dG%len_lat = oG%len_lat ; dG%len_lon = oG%len_lon

src/framework/MOM_dyn_horgrid.F90

@@ -181,6 +181,10 @@ module MOM_dyn_horgrid
 
   ! These parameters are run-time parameters that are used during some
   ! initialization routines (but not all)
+  real :: grid_unit_to_L !< A factor that converts a the geoLat and geoLon variables and related
+                        !! variables like len_lat and len_lon into rescaled horizontal distance
+                        !! units on a Cartesian grid, in [L km ~> 1000] or [L m-1 ~> 1] or
+                        !! is 0 for a non-Cartesian grid.
   real :: south_lat     !< The latitude (or y-coordinate) of the first v-line [degrees_N] or [km] or [m]
   real :: west_lon      !< The longitude (or x-coordinate) of the first u-line [degrees_E] or [km] or [m]
   real :: len_lat       !< The latitudinal (or y-coord) extent of physical domain [degrees_N] or [km] or [m]
@@ -400,6 +404,7 @@ subroutine rotate_dyn_horgrid(G_in, G, US, turns)
   G%len_lon = G_in%len_lon
 
   ! Rotation-invariant fields
+  G%grid_unit_to_L = G_in%grid_unit_to_L
   G%areaT_global = G_in%areaT_global
   G%IareaT_global = G_in%IareaT_global
   G%Rad_Earth = G_in%Rad_Earth

src/ice_shelf/MOM_ice_shelf_diag_mediator.F90

@@ -132,6 +132,7 @@ subroutine set_IS_axes_info(G, param_file, diag_cs, axes_set_name)
 
   G%x_axis_units = "degrees_E" ; G%y_axis_units = "degrees_N"
   G%x_ax_unit_short = "degrees_E" ; G%y_ax_unit_short = "degrees_N"
+  G%grid_unit_to_L = 0.0
 
   if (index(lowercase(trim(grid_config)),"cartesian") > 0) then
     ! This is a cartesian grid, and may have different axis units.
@@ -145,9 +146,11 @@ subroutine set_IS_axes_info(G, param_file, diag_cs, axes_set_name)
     if (units_temp(1:1) == 'k') then
       G%x_axis_units = "kilometers" ; G%y_axis_units = "kilometers"
       G%x_ax_unit_short = "km" ; G%y_ax_unit_short = "km"
+      G%grid_unit_to_L = 1000.0*diag_cs%US%m_to_L
     elseif (units_temp(1:1) == 'm') then
       G%x_axis_units = "meters" ; G%y_axis_units = "meters"
       G%x_ax_unit_short = "m" ; G%y_ax_unit_short = "m"
+      G%grid_unit_to_L = diag_cs%US%m_to_L
     endif
     call log_param(param_file, mdl, "explicit AXIS_UNITS", G%x_axis_units)
   else

src/initialization/MOM_grid_initialize.F90

@@ -84,7 +84,7 @@ subroutine set_grid_metrics(G, param_file, US)
   ! These are defaults that may be changed in the next select block.
   G%x_axis_units = "degrees_east" ; G%y_axis_units = "degrees_north"
   G%x_ax_unit_short = "degrees_E" ; G%y_ax_unit_short = "degrees_N"
-
+  G%grid_unit_to_L = 0.0
   G%Rad_Earth_L = -1.0*US%m_to_L ; G%len_lat = 0.0 ; G%len_lon = 0.0
   select case (trim(config))
     case ("mosaic");    call set_grid_metrics_from_mosaic(G, param_file, US)
@@ -251,6 +251,11 @@ subroutine set_grid_metrics_from_mosaic(G, param_file, US)
     G%geoLatCv(i,J) = tmpZ(i2-1,j2)
   enddo ; enddo
 
+  ! This routine could be modified to support the use of a mosaic using Cartesian grid coordinates,
+  ! in which case the values of G%x_axis_units, G%y_axis_units and G%grid_unit_to_L would need to be
+  ! reset appropriately here, but this option has not yet been implemented, and the grid coordinates
+  ! are assumed to be degrees of longitude and latitude.
+
   ! Read DX,DY from the supergrid
   tmpU(:,:) = 0. ; tmpV(:,:) = 0.
   call MOM_read_data(filename, 'dx', tmpV, SGdom, position=NORTH_FACE, scale=US%m_to_L)
@@ -440,9 +445,11 @@ subroutine set_grid_metrics_cartesian(G, param_file, US)
   enddo
 
   if (units_temp(1:1) == 'k') then ! Axes are measured in km.
+    G%grid_unit_to_L = 1000.0*US%m_to_L
     dx_everywhere = 1000.0*US%m_to_L * G%len_lon / (REAL(niglobal))
     dy_everywhere = 1000.0*US%m_to_L * G%len_lat / (REAL(njglobal))
   elseif (units_temp(1:1) == 'm') then ! Axes are measured in m.
+    G%grid_unit_to_L = US%m_to_L
     dx_everywhere = US%m_to_L*G%len_lon / (REAL(niglobal))
     dy_everywhere = US%m_to_L*G%len_lat / (REAL(njglobal))
   else ! Axes are measured in degrees of latitude and longitude.

src/initialization/MOM_shared_initialization.F90

@@ -485,7 +485,6 @@ subroutine set_rotation_beta_plane(f, G, param_file, US)
   real    :: f_0    ! The reference value of the Coriolis parameter [T-1 ~> s-1]
   real    :: beta   ! The meridional gradient of the Coriolis parameter [T-1 L-1 ~> s-1 m-1]
   real    :: beta_lat_ref ! The reference latitude for the beta plane [degrees_N] or [km] or [m]
-  real    :: Rad_Earth_L  ! The radius of the planet in rescaled units [L ~> m]
   real    :: y_scl  ! A scaling factor from the units of latitude [L lat-1 ~> m lat-1]
   real    :: PI     ! The ratio of the circumference of a circle to its diameter [nondim]
   character(len=40)  :: mdl = "set_rotation_beta_plane" ! This subroutine's name.
@@ -503,18 +502,16 @@ subroutine set_rotation_beta_plane(f, G, param_file, US)
   call get_param(param_file, mdl, "AXIS_UNITS", axis_units, default="degrees")
 
   PI = 4.0*atan(1.0)
+  y_scl = G%grid_unit_to_L
+  if (G%grid_unit_to_L <= 0.0) y_scl = PI * G%Rad_Earth_L / 180.
+
   select case (axis_units(1:1))
     case ("d")
-      call get_param(param_file, mdl, "RAD_EARTH", Rad_Earth_L, &
-                   "The radius of the Earth.", units="m", default=6.378e6, scale=US%m_to_L)
       beta_lat_ref_units = "degrees"
-      y_scl = PI * Rad_Earth_L / 180.
     case ("k")
       beta_lat_ref_units = "kilometers"
-      y_scl = 1.0e3 * US%m_to_L
     case ("m")
       beta_lat_ref_units = "meters"
-      y_scl = 1.0 * US%m_to_L
     case default ; call MOM_error(FATAL, &
       " set_rotation_beta_plane: unknown AXIS_UNITS = "//trim(axis_units))
   end select

src/initialization/MOM_state_initialization.F90

@@ -1564,7 +1564,10 @@ subroutine initialize_velocity_circular(u, v, G, GV, US, param_file, just_read)
 
   if (just_read) return ! All run-time parameters have been read, so return.
 
-  dpi=acos(0.0)*2.0 ! pi
+  if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "MOM_state_initialization.F90: "//&
+          "initialize_velocity_circular() is only set to work with Cartesian axis units.")
+
+  dpi = acos(0.0)*2.0 ! pi
 
   do k=1,nz ; do j=js,je ; do I=Isq,Ieq
     psi1 = my_psi(I,j)
@@ -1591,7 +1594,7 @@ real function my_psi(ig,jg)
     r = sqrt( (x**2) + (y**2) ) ! Circular stream function is a function of radius only
     r = min(1.0, r) ! Flatten stream function in corners of box
     my_psi = 0.5*(1.0 - cos(dpi*r))
-    my_psi = my_psi * (circular_max_u * G%US%m_to_L*G%len_lon*1e3 / dpi) ! len_lon is in km
+    my_psi = my_psi * (circular_max_u * G%len_lon * G%grid_unit_to_L / dpi) ! len_lon is in km
   end function my_psi
 
 end subroutine initialize_velocity_circular

src/user/DOME_initialization.F90

@@ -49,18 +49,28 @@ subroutine DOME_initialize_topography(D, G, param_file, max_depth, US)
   real :: min_depth  ! The minimum ocean depth [Z ~> m]
   real :: shelf_depth ! The ocean depth on the shelf in the DOME configuration [Z ~> m]
   real :: slope      ! The bottom slope in the DOME configuration [Z L-1 ~> nondim]
-  real :: shelf_edge_lat ! The latitude of the edge of the topographic shelf [km]
-  real :: inflow_lon ! The edge longitude of the DOME inflow [km]
-  real :: inflow_width ! The longitudinal width of the DOME inflow channel [km]
-  real :: km_to_L    ! The conversion factor from the units of latitude to L [L km-1 ~> 1e3]
+  real :: shelf_edge_lat ! The latitude of the edge of the topographic shelf in the same units as geolat, often [km]
+  real :: inflow_lon ! The edge longitude of the DOME inflow in the same units as geolon, often [km]
+  real :: inflow_width ! The longitudinal width of the DOME inflow channel in the same units as geolat, often [km]
+  real :: km_to_grid_unit ! The conversion factor from km to the units of latitude often 1 [nondim],
+                          ! but this could be 1000 [m km-1]
   ! This include declares and sets the variable "version".
 # include "version_variable.h"
   character(len=40)  :: mdl = "DOME_initialize_topography" ! This subroutine's name.
   integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
 
-  km_to_L = 1.0e3*US%m_to_L
+  if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "DOME_initialization: "//&
+          "DOME_initialize_topography is only set to work with Cartesian axis units.")
+  if (abs(G%grid_unit_to_L*US%L_to_m - 1000.0) < 1.0e-3) then ! The grid latitudes are in km.
+    km_to_grid_unit = 1.0
+  elseif (abs(G%grid_unit_to_L*US%L_to_m - 1.0) < 1.0e-6) then ! The grid latitudes are in m.
+    km_to_grid_unit = 1000.0
+  else
+    call MOM_error(FATAL, "DOME_initialization: "//&
+          "DOME_initialize_topography is not recognizing the value of G%grid_unit_to_L.")
+  endif
 
   call MOM_mesg("  DOME_initialization.F90, DOME_initialize_topography: setting topography", 5)
 
@@ -75,15 +85,16 @@ subroutine DOME_initialize_topography(D, G, param_file, max_depth, US)
                  default=600.0, units="m", scale=US%m_to_Z)
   call get_param(param_file, mdl, "DOME_SHELF_EDGE_LAT", shelf_edge_lat, &
                  "The latitude of the shelf edge in the DOME configuration.", &
-                 default=600.0, units="km")
+                 default=600.0, units="km", scale=km_to_grid_unit)
   call get_param(param_file, mdl, "DOME_INFLOW_LON", inflow_lon, &
-                 "The edge longitude of the DOME inflow.", units="km", default=1000.0)
+                 "The edge longitude of the DOME inflow.", units="km", default=1000.0, scale=km_to_grid_unit)
   call get_param(param_file, mdl, "DOME_INFLOW_WIDTH", inflow_width, &
-                 "The longitudinal width of the DOME inflow channel.", units="km", default=100.0)
+                 "The longitudinal width of the DOME inflow channel.", &
+                 units="km", default=100.0, scale=km_to_grid_unit)
 
   do j=js,je ; do i=is,ie
     if (G%geoLatT(i,j) < shelf_edge_lat) then
-      D(i,j) = min(shelf_depth - slope * (G%geoLatT(i,j)-shelf_edge_lat)*km_to_L, max_depth)
+      D(i,j) = min(shelf_depth - slope * (G%geoLatT(i,j)-shelf_edge_lat)*G%grid_unit_to_L, max_depth)
     else
       if ((G%geoLonT(i,j) > inflow_lon) .AND. (G%geoLonT(i,j) < inflow_lon+inflow_width)) then
         D(i,j) = shelf_depth
@@ -177,7 +188,6 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp)
   real :: min_depth    ! The minimum depth at which to apply damping [Z ~> m]
   real :: damp_W, damp_E ! Damping rates in the western and eastern sponges [T-1 ~> s-1]
   real :: peak_damping ! The maximum sponge damping rates as the edges [T-1 ~> s-1]
-  real :: km_to_L      ! The conversion factor from the units of longitude to L [L km-1 ~> 1e3]
   real :: edge_dist    ! The distance to an edge [L ~> m]
   real :: sponge_width ! The width of the sponges [L ~> m]
   character(len=40)  :: mdl = "DOME_initialize_sponges" ! This subroutine's name.
@@ -186,7 +196,8 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp)
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
 
-  km_to_L = 1.0e3*US%m_to_L
+  if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "DOME_initialization: "//&
+          "DOME_initialize_sponges is only set to work with Cartesian axis units.")
 
   !   Set up sponges for the DOME configuration
   call get_param(PF, mdl, "MINIMUM_DEPTH", min_depth, &
@@ -196,15 +207,15 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp)
                  default=10.0, units="day-1", scale=1.0/(86400.0*US%s_to_T))
   call get_param(PF, mdl, "DOME_SPONGE_WIDTH", sponge_width, &
                  "The width of the the DOME sponges.", &
-                 default=200.0, units="km", scale=km_to_L)
+                 default=200.0, units="km", scale=1.0e3*US%m_to_L)
 
   ! Here the inverse damping time [T-1 ~> s-1], is set. Set Idamp to 0 wherever
   ! there is no sponge, and the subroutines that are called will automatically
   ! set up the sponges only where Idamp is positive and mask2dT is 1.
 
   Idamp(:,:) = 0.0
   do j=js,je ; do i=is,ie ; if (depth_tot(i,j) > min_depth) then
-    edge_dist = (G%geoLonT(i,j) - G%west_lon) * km_to_L
+    edge_dist = (G%geoLonT(i,j) - G%west_lon) * G%grid_unit_to_L
     if (edge_dist < 0.5*sponge_width) then
       damp_W = peak_damping
     elseif (edge_dist < sponge_width) then
@@ -213,7 +224,7 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp)
       damp_W = 0.0
     endif
 
-    edge_dist = ((G%len_lon + G%west_lon) - G%geoLonT(i,j)) * km_to_L
+    edge_dist = ((G%len_lon + G%west_lon) - G%geoLonT(i,j)) * G%grid_unit_to_L
     if (edge_dist < 0.5*sponge_width) then
       damp_E = peak_damping
     elseif (edge_dist < sponge_width) then
@@ -328,10 +339,12 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg)
                             ! properties [T-1 ~> s-1]
   real :: g_prime_tot       ! The reduced gravity across all layers [L2 Z-1 T-2 ~> m s-2]
   real :: Def_Rad           ! The deformation radius, based on fluid of thickness D_edge [L ~> m]
-  real :: inflow_lon        ! The edge longitude of the DOME inflow [km]
+  real :: inflow_lon        ! The edge longitude of the DOME inflow in the same units as geolon, often [km]
   real :: I_Def_Rad         ! The inverse of the deformation radius in the same units as G%geoLon [km-1]
   real :: Ri_trans          ! The shear Richardson number in the transition
                             ! region of the specified shear profile [nondim]
+  real :: km_to_grid_unit   ! The conversion factor from km to the units of latitude often 1 [nondim],
+                            ! but this could be 1000 [m km-1]
   character(len=32)  :: name ! The name of a tracer field.
   character(len=40)  :: mdl = "DOME_set_OBC_data" ! This subroutine's name.
   integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, m, nz, ntherm, ntr_id
@@ -343,6 +356,17 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg)
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
   IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
 
+  if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "DOME_initialization: "//&
+          "DOME_initialize_topography is only set to work with Cartesian axis units.")
+  if (abs(G%grid_unit_to_L*US%L_to_m - 1000.0) < 1.0e-3) then ! The grid latitudes are in km.
+    km_to_grid_unit = 1.0
+  elseif (abs(G%grid_unit_to_L*US%L_to_m - 1.0) < 1.0e-6) then ! The grid latitudes are in m.
+    km_to_grid_unit = 1000.0
+  else
+    call MOM_error(FATAL, "DOME_initialization: "//&
+          "DOME_initialize_topography is not recognizing the value of G%grid_unit_to_L.")
+  endif
+
   call get_param(PF, mdl, "DOME_INFLOW_THICKNESS", D_edge, &
                  "The thickness of the dense DOME inflow at the inner edge.", &
                  default=300.0, units="m", scale=US%m_to_Z)
@@ -362,7 +386,7 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg)
                  "The value of the Coriolis parameter that is used to determine the DOME "//&
                  "inflow properties.", units="s-1", default=f_0*US%s_to_T, scale=US%T_to_s)
   call get_param(PF, mdl, "DOME_INFLOW_LON", inflow_lon, &
-                 "The edge longitude of the DOME inflow.", units="km", default=1000.0)
+                 "The edge longitude of the DOME inflow.", units="km", default=1000.0, scale=km_to_grid_unit)
   if (associated(tv%S) .or. associated(tv%T)) then
     call get_param(PF, mdl, "S_REF", S_ref, &
                  units="ppt", default=35.0, scale=US%ppt_to_S, do_not_log=.true.)
@@ -383,7 +407,9 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg)
     tr_0 = (-D_edge*sqrt(D_edge*g_prime_tot)*0.5*Def_Rad) * (Rlay_Ref + 0.5*Rlay_range) * GV%RZ_to_H
   endif
 
-  I_Def_Rad = 1.0 / (1.0e-3*US%L_to_m*Def_Rad)
+  I_Def_Rad = 1.0 / ((1.0e-3*US%L_to_m*km_to_grid_unit) * Def_Rad)
+  ! This is mathematically equivalent to
+  ! I_Def_Rad = G%grid_unit_to_L / Def_Rad
 
   if (OBC%number_of_segments /= 1) then
     call MOM_error(WARNING, 'Error in DOME OBC segment setup', .true.)

src/user/ISOMIP_initialization.F90

@@ -59,7 +59,6 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US)
   real :: ly              ! domain width (across ice flow) [L ~> m]
   real :: bx, by          ! The x- and y- contributions to the bathymetric profiles at a point [Z ~> m]
   real :: xtil            ! x-positon normalized by the characteristic along-flow length scale [nondim]
-  real :: km_to_L         ! The conversion factor from the axis units to L [L km-1 ~> 1e3]
   logical :: is_2D        ! If true, use a 2D setup
   ! This include declares and sets the variable "version".
 # include "version_variable.h"
@@ -93,16 +92,17 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US)
                  "Characteristic width of the side walls of the channel in the ISOMIP configuration.", &
                  units="m", default=4.0e3, scale=US%m_to_L)
 
-  km_to_L = 1.0e3*US%m_to_L
+  if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "ISOMIP_initialization.F90: " //&
+          "ISOMIP_initialize_topography is only set to work with Cartesian axis units.")
 
   ! The following variables should be transformed into runtime parameters.
   b0 = -150.0*US%m_to_Z ; b2 = -728.8*US%m_to_Z ; b4 = 343.91*US%m_to_Z ; b6 = -50.57*US%m_to_Z
 
   if (is_2D) then
     do j=js,je ; do i=is,ie
       ! For the 2D setup take a slice through the middle of the domain
-      xtil = G%geoLonT(i,j)*km_to_L / xbar
-      !xtil = 450.*km_to_L / xbar
+      xtil = G%geoLonT(i,j)*G%grid_unit_to_L / xbar
+      ! xtil = 450.0e3*US%m_to_L / xbar
       bx = b0 + b2*xtil**2 + b4*xtil**4 + b6*xtil**6
 
       by = 2.0 * dc / (1.0 + exp(2.0*wc / fc))
@@ -117,17 +117,17 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US)
       ! 3D setup
       ! ===== TEST =====
       !if (G%geoLonT(i,j)<500.) then
-      !  xtil = 500.*km_to_L / xbar
+      !  xtil = 500.0e3*US%m_to_L / xbar
       !else
-      !  xtil = G%geoLonT(i,j)*km_to_L / xbar
+      !  xtil = G%geoLonT(i,j)*G%grid_unit_to_L / xbar
       !endif
       ! ===== TEST =====
 
-      xtil = G%geoLonT(i,j)*km_to_L / xbar
+      xtil = G%geoLonT(i,j)*G%grid_unit_to_L / xbar
 
       bx = b0 + b2*xtil**2 + b4*xtil**4 + b6*xtil**6
-      by = (dc / (1.0 + exp(-2.*(G%geoLatT(i,j)*km_to_L - 0.5*ly - wc) / fc))) + &
-           (dc / (1.0 + exp(2.*(G%geoLatT(i,j)*km_to_L - 0.5*ly + wc) / fc)))
+      by = (dc / (1.0 + exp(-2.*(G%geoLatT(i,j)*G%grid_unit_to_L - 0.5*ly - wc) / fc))) + &
+           (dc / (1.0 + exp(2.*(G%geoLatT(i,j)*G%grid_unit_to_L - 0.5*ly + wc) / fc)))
 
       D(i,j) = -max(bx+by, -bmax)
       if (D(i,j) > max_depth) D(i,j) = max_depth