MOM_shared_initialization.F90

!> Code that initializes fixed aspects of the model grid, such as horizontal
!! grid metrics, topography and Coriolis, and can be shared between components.
module mom_shared_initialization

! This file is part of MOM6. See LICENSE.md for the license.

use mom_coms, only : max_across_pes, reproducing_sum
use mom_domains, only : pass_var, pass_vector, sum_across_pes, broadcast
use mom_domains, only : root_pe, to_all, scalar_pair, cgrid_ne, agrid
use mom_dyn_horgrid, only : dyn_horgrid_type
use mom_error_handler, only : mom_mesg, mom_error, fatal, warning, is_root_pe
use mom_error_handler, only : calltree_enter, calltree_leave, calltree_waypoint
use mom_file_parser, only : get_param, log_param, param_file_type, log_version
use mom_io, only : close_file, create_file, fieldtype, file_exists
use mom_io, only : mom_read_data, mom_read_vector, single_file, multiple
use mom_io, only : slasher, vardesc, write_field, var_desc
use mom_string_functions, only : uppercase
use mom_unit_scaling, only : unit_scale_type

use netcdf

implicit none ; private

public mom_shared_init_init
public mom_initialize_rotation, mom_calculate_grad_coriolis
public initialize_topography_from_file, apply_topography_edits_from_file
public initialize_topography_named, limit_topography, diagnosemaximumdepth
public set_rotation_planetary, set_rotation_beta_plane, initialize_grid_rotation_angle
public reset_face_lengths_named, reset_face_lengths_file, reset_face_lengths_list
public read_face_length_list, set_velocity_depth_max, set_velocity_depth_min
public compute_global_grid_integrals, write_ocean_geometry_file

! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.

contains

! -----------------------------------------------------------------------------
!> MOM_shared_init_init just writes the code version.
subroutine mom_shared_init_init(PF)
  type(param_file_type),   intent(in)    :: pf   !< A structure indicating the open file
                                                 !! to parse for model parameter values.

  character(len=40)  :: mdl = "MOM_shared_initialization" ! This module's name.

! This include declares and sets the variable "version".
#include "version_variable.h"
  call log_version(pf, mdl, version, &
   "Sharable code to initialize time-invariant fields, like bathymetry and Coriolis parameters.")

end subroutine mom_shared_init_init
! -----------------------------------------------------------------------------

!> MOM_initialize_rotation makes the appropriate call to set up the Coriolis parameter.
subroutine mom_initialize_rotation(f, G, PF, US)
  type(dyn_horgrid_type),                       intent(in)  :: g  !< The dynamic horizontal grid type
  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), intent(out) :: f  !< The Coriolis parameter [T-1 ~> s-1]
  type(param_file_type),                        intent(in)  :: pf !< Parameter file structure
  type(unit_scale_type),              optional, intent(in)  :: us !< A dimensional unit scaling type

!   This subroutine makes the appropriate call to set up the Coriolis parameter.
! This is a separate subroutine so that it can be made public and shared with
! the ice-sheet code or other components.
! Set up the Coriolis parameter, f, either analytically or from file.
  character(len=40)  :: mdl = "MOM_initialize_rotation" ! This subroutine's name.
  character(len=200) :: config

  call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
  call get_param(pf, mdl, "ROTATION", config, &
                 "This specifies how the Coriolis parameter is specified: \n"//&
                 " \t 2omegasinlat - Use twice the planetary rotation rate \n"//&
                 " \t\t times the sine of latitude.\n"//&
                 " \t betaplane - Use a beta-plane or f-plane.\n"//&
                 " \t USER - call a user modified routine.", &
                 default="2omegasinlat")
  select case (trim(config))
    case ("2omegasinlat"); call set_rotation_planetary(f, g, pf, us)
    case ("beta"); call set_rotation_beta_plane(f, g, pf, us)
    case ("betaplane"); call set_rotation_beta_plane(f, g, pf, us)
   !case ("nonrotating") ! Note from AJA: Missing case?
    case default ; call mom_error(fatal,"MOM_initialize: "// &
      "Unrecognized rotation setup "//trim(config))
  end select
  call calltree_leave(trim(mdl)//'()')
end subroutine mom_initialize_rotation

!> Calculates the components of grad f (Coriolis parameter)
subroutine mom_calculate_grad_coriolis(dF_dx, dF_dy, G, US)
  type(dyn_horgrid_type),             intent(inout) :: g !< The dynamic horizontal grid type
  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
                                      intent(out)   :: df_dx !< x-component of grad f [T-1 L-1 ~> s-1 m-1]
  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
                                      intent(out)   :: df_dy !< y-component of grad f [T-1 L-1 ~> s-1 m-1]
  type(unit_scale_type),    optional, intent(in)    :: us !< A dimensional unit scaling type
  ! Local variables
  integer :: i,j
  real :: m_to_l  ! A unit conversion factor [L m-1 ~> nondim]
  real :: f1, f2

  m_to_l = 1.0 ; if (present(us)) m_to_l = us%m_to_L

  if ((lbound(g%CoriolisBu,1) > g%isc-1) .or. &
      (lbound(g%CoriolisBu,2) > g%isc-1)) then
    ! The gradient of the Coriolis parameter can not be calculated with this grid.
    df_dx(:,:) = 0.0 ; df_dy(:,:) = 0.0
    return
  endif

  do j=g%jsc, g%jec ; do i=g%isc, g%iec
    f1 = 0.5*( g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1) )
    f2 = 0.5*( g%CoriolisBu(i-1,j) + g%CoriolisBu(i-1,j-1) )
    df_dx(i,j) = g%IdxT(i,j) * ( f1 - f2 )
    f1 = 0.5*( g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j) )
    f2 = 0.5*( g%CoriolisBu(i,j-1) + g%CoriolisBu(i-1,j-1) )
    df_dy(i,j) = g%IdyT(i,j) * ( f1 - f2 )
  enddo ; enddo
  call pass_vector(df_dx, df_dy, g%Domain, stagger=agrid)

end subroutine mom_calculate_grad_coriolis

!> Return the global maximum ocean bottom depth in the same units as the input depth.
function diagnosemaximumdepth(D, G)
  type(dyn_horgrid_type),  intent(in) :: g !< The dynamic horizontal grid type
  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
                           intent(in) :: d !< Ocean bottom depth in m or Z
  real :: diagnosemaximumdepth             !< The global maximum ocean bottom depth in m or Z
  ! Local variables
  integer :: i,j
  diagnosemaximumdepth = d(g%isc,g%jsc)
  do j=g%jsc, g%jec ; do i=g%isc, g%iec
    diagnosemaximumdepth = max(diagnosemaximumdepth,d(i,j))
  enddo ; enddo
  call max_across_pes(diagnosemaximumdepth)
end function diagnosemaximumdepth


!> Read gridded depths from file
subroutine initialize_topography_from_file(D, G, param_file, US)
  type(dyn_horgrid_type),           intent(in)  :: g !< The dynamic horizontal grid type
  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
                                    intent(out) :: d !< Ocean bottom depth in m or Z if US is present
  type(param_file_type),            intent(in)  :: param_file !< Parameter file structure
  type(unit_scale_type),  optional, intent(in)  :: us !< A dimensional unit scaling type
  ! Local variables
  real :: m_to_z  ! A dimensional rescaling factor.
  character(len=200) :: filename, topo_file, inputdir ! Strings for file/path
  character(len=200) :: topo_varname                  ! Variable name in file
  character(len=40)  :: mdl = "initialize_topography_from_file" ! This subroutine's name.

  call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")

  m_to_z = 1.0 ; if (present(us)) m_to_z = us%m_to_Z

  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  call get_param(param_file, mdl, "TOPO_FILE", topo_file, &
                 "The file from which the bathymetry is read.", &
                 default="topog.nc")
  call get_param(param_file, mdl, "TOPO_VARNAME", topo_varname, &
                 "The name of the bathymetry variable in TOPO_FILE.", &
                 default="depth")

  filename = trim(inputdir)//trim(topo_file)
  call log_param(param_file, mdl, "INPUTDIR/TOPO_FILE", filename)

  if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
       " initialize_topography_from_file: Unable to open "//trim(filename))

  d(:,:) = -9.e30*m_to_z ! Initializing to a very large negative depth (tall mountains) everywhere
                         ! before reading from a file should do nothing. However, in the instance of
                         ! masked-out PEs, halo regions are not updated when a processor does not
                         ! exist. We need to ensure the depth in masked-out PEs appears to be that
                         ! of land so this line does that in the halo regions. For non-masked PEs
                         ! the halo region is filled properly with a later pass_var().
  call mom_read_data(filename, trim(topo_varname), d, g%Domain, scale=m_to_z)

  call apply_topography_edits_from_file(d, g, param_file, us)

  call calltree_leave(trim(mdl)//'()')
end subroutine initialize_topography_from_file

!> Applies a list of topography overrides read from a netcdf file
subroutine apply_topography_edits_from_file(D, G, param_file, US)
  type(dyn_horgrid_type),           intent(in)    :: g !< The dynamic horizontal grid type
  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
                                    intent(inout) :: d !< Ocean bottom depth in m or Z if US is present
  type(param_file_type),            intent(in)    :: param_file !< Parameter file structure
  type(unit_scale_type),  optional, intent(in)    :: us !< A dimensional unit scaling type

  ! Local variables
  real :: m_to_z  ! A dimensional rescaling factor.
  character(len=200) :: topo_edits_file, inputdir ! Strings for file/path
  character(len=40)  :: mdl = "apply_topography_edits_from_file" ! This subroutine's name.
  integer :: n_edits, n, ashape(5), i, j, ncid, id, ncstatus, iid, jid, zid
  integer, dimension(:), allocatable :: ig, jg
  real, dimension(:), allocatable :: new_depth

  call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")

  m_to_z = 1.0 ; if (present(us)) m_to_z = us%m_to_Z

  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  call get_param(param_file, mdl, "TOPO_EDITS_FILE", topo_edits_file, &
                 "The file from which to read a list of i,j,z topography overrides.", &
                 default="")

  if (len_trim(topo_edits_file)==0) return

  topo_edits_file = trim(inputdir)//trim(topo_edits_file)
  if (.not.file_exists(topo_edits_file, g%Domain)) call mom_error(fatal, &
     'initialize_topography_from_file: Unable to open '//trim(topo_edits_file))

  ncstatus = nf90_open(trim(topo_edits_file), nf90_nowrite, ncid)
  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
                                'Failed to open '//trim(topo_edits_file))

  ! Get nEdits
  ncstatus = nf90_inq_dimid(ncid, 'nEdits', id)
  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
                                'Failed to inq_dimid nEdits for '//trim(topo_edits_file))
  ncstatus = nf90_inquire_dimension(ncid, id, len=n_edits)
  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
                                'Failed to inquire_dimension nEdits for '//trim(topo_edits_file))

  ! Read ni
  ncstatus = nf90_inq_varid(ncid, 'ni', id)
  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
                                'Failed to inq_varid ni for '//trim(topo_edits_file))
  ncstatus = nf90_get_var(ncid, id, i)
  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
                              'Failed to get_var ni for '//trim(topo_edits_file))
  if (i /= g%ieg) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
                              'Incompatible i-dimension of grid in '//trim(topo_edits_file))

  ! Read nj
  ncstatus = nf90_inq_varid(ncid, 'nj', id)
  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
                                'Failed to inq_varid nj for '//trim(topo_edits_file))
  ncstatus = nf90_get_var(ncid, id, j)
  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
                              'Failed to get_var nj for '//trim(topo_edits_file))
  if (j /= g%jeg) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
                              'Incompatible j-dimension of grid in '//trim(topo_edits_file))

  ! Read iEdit
  ncstatus = nf90_inq_varid(ncid, 'iEdit', id)
  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
                                'Failed to inq_varid iEdit for '//trim(topo_edits_file))
  allocate(ig(n_edits))
  ncstatus = nf90_get_var(ncid, id, ig)
  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
                              'Failed to get_var iEdit for '//trim(topo_edits_file))

  ! Read jEdit
  ncstatus = nf90_inq_varid(ncid, 'jEdit', id)
  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
                                'Failed to inq_varid jEdit for '//trim(topo_edits_file))
  allocate(jg(n_edits))
  ncstatus = nf90_get_var(ncid, id, jg)
  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
                              'Failed to get_var jEdit for '//trim(topo_edits_file))

  ! Read zEdit
  ncstatus = nf90_inq_varid(ncid, 'zEdit', id)
  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
                                'Failed to inq_varid zEdit for '//trim(topo_edits_file))
  allocate(new_depth(n_edits))
  ncstatus = nf90_get_var(ncid, id, new_depth)
  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
                              'Failed to get_var zEdit for '//trim(topo_edits_file))

  ! Close file
  ncstatus = nf90_close(ncid)
  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
                                'Failed to close '//trim(topo_edits_file))

  do n = 1, n_edits
    i = ig(n) - g%isd_global + 2 ! +1 for python indexing and +1 for ig-isd_global+1
    j = jg(n) - g%jsd_global + 2
    if (i>=g%isc .and. i<=g%iec .and. j>=g%jsc .and. j<=g%jec) then
      if (new_depth(n)/=0.) then
        write(*,'(a,3i5,f8.2,a,f8.2,2i4)') &
          'Ocean topography edit: ',n,ig(n),jg(n),d(i,j)/m_to_z,'->',abs(new_depth(n)),i,j
        d(i,j) = abs(m_to_z*new_depth(n)) ! Allows for height-file edits (i.e. converts negatives)
      else
        call mom_error(fatal, ' apply_topography_edits_from_file: '//&
          "A zero depth edit would change the land mask and is not allowed in"//trim(topo_edits_file))
      endif
    endif
  enddo

  deallocate( ig, jg, new_depth )

  call calltree_leave(trim(mdl)//'()')
end subroutine apply_topography_edits_from_file

!> initialize the bathymetry based on one of several named idealized configurations
subroutine initialize_topography_named(D, G, param_file, topog_config, max_depth, US)
  type(dyn_horgrid_type),           intent(in)  :: g !< The dynamic horizontal grid type
  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
                                    intent(out) :: d !< Ocean bottom depth in m or Z if US is present
  type(param_file_type),            intent(in)  :: param_file !< Parameter file structure
  character(len=*),                 intent(in)  :: topog_config !< The name of an idealized
                                                              !! topographic configuration
  real,                             intent(in)  :: max_depth  !< Maximum depth of model in the units of D
  type(unit_scale_type),  optional, intent(in)  :: us !< A dimensional unit scaling type

  ! This subroutine places the bottom depth in m into D(:,:), shaped according to the named config.

  ! Local variables
  real :: m_to_z               ! A dimensional rescaling factor.
  real :: min_depth            ! The minimum depth [Z ~> m].
  real :: pi                   ! 3.1415926... calculated as 4*atan(1)
  real :: d0                   ! A constant to make the maximum  basin depth MAXIMUM_DEPTH.
  real :: expdecay             ! A decay scale of associated with the sloping boundaries [m].
  real :: dedge                ! The depth [Z ~> m], at the basin edge
! real :: south_lat, west_lon, len_lon, len_lat, Rad_earth
  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
  character(len=40)  :: mdl = "initialize_topography_named" ! This subroutine's name.
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

  call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
  call mom_mesg("  MOM_shared_initialization.F90, initialize_topography_named: "//&
                 "TOPO_CONFIG = "//trim(topog_config), 5)

  m_to_z = 1.0 ; if (present(us)) m_to_z = us%m_to_Z

  call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
                 "The minimum depth of the ocean.", units="m", default=0.0, scale=m_to_z)
  if (max_depth<=0.) call mom_error(fatal,"initialize_topography_named: "// &
      "MAXIMUM_DEPTH has a non-sensical value! Was it set?")

  if (trim(topog_config) /= "flat") then
    call get_param(param_file, mdl, "EDGE_DEPTH", dedge, &
                   "The depth at the edge of one of the named topographies.", &
                   units="m", default=100.0, scale=m_to_z)
!   call get_param(param_file, mdl, "SOUTHLAT", south_lat, &
!                  "The southern latitude of the domain.", units="degrees", &
!                  fail_if_missing=.true.)
!   call get_param(param_file, mdl, "LENLAT", len_lat, &
!                  "The latitudinal length of the domain.", units="degrees", &
!                  fail_if_missing=.true.)
!   call get_param(param_file, mdl, "WESTLON", west_lon, &
!                  "The western longitude of the domain.", units="degrees", &
!                  default=0.0)
!   call get_param(param_file, mdl, "LENLON", len_lon, &
!                  "The longitudinal length of the domain.", units="degrees", &
!                  fail_if_missing=.true.)
!   call get_param(param_file, mdl, "RAD_EARTH", Rad_Earth, &
!                  "The radius of the Earth.", units="m", default=6.378e6)
    call get_param(param_file, mdl, "TOPOG_SLOPE_SCALE", expdecay, &
                   "The exponential decay scale used in defining some of "//&
                   "the named topographies.", units="m", default=400000.0)
  endif


  pi = 4.0*atan(1.0)

  if (trim(topog_config) == "flat") then
    do i=is,ie ; do j=js,je ; d(i,j) = max_depth ; enddo ; enddo
  elseif (trim(topog_config) == "spoon") then
    d0 = (max_depth - dedge) / &
             ((1.0 - exp(-0.5*g%len_lat*g%Rad_earth*pi/(180.0 *expdecay))) * &
              (1.0 - exp(-0.5*g%len_lat*g%Rad_earth*pi/(180.0 *expdecay))))
    do i=is,ie ; do j=js,je
  !  This sets a bowl shaped (sort of) bottom topography, with a       !
  !  maximum depth of max_depth.                                   !
      d(i,j) =  dedge + d0 * &
             (sin(pi * (g%geoLonT(i,j) - (g%west_lon)) / g%len_lon) * &
           (1.0 - exp((g%geoLatT(i,j) - (g%south_lat+g%len_lat))*g%Rad_earth*pi / &
                      (180.0*expdecay)) ))
    enddo ; enddo
  elseif (trim(topog_config) == "bowl") then
    d0 = (max_depth - dedge) / &
             ((1.0 - exp(-0.5*g%len_lat*g%Rad_earth*pi/(180.0 *expdecay))) * &
              (1.0 - exp(-0.5*g%len_lat*g%Rad_earth*pi/(180.0 *expdecay))))

  !  This sets a bowl shaped (sort of) bottom topography, with a
  !  maximum depth of max_depth.
    do i=is,ie ; do j=js,je
      d(i,j) =  dedge + d0 * &
             (sin(pi * (g%geoLonT(i,j) - g%west_lon) / g%len_lon) * &
             ((1.0 - exp(-(g%geoLatT(i,j) - g%south_lat)*g%Rad_Earth*pi/ &
                          (180.0*expdecay))) * &
             (1.0 - exp((g%geoLatT(i,j) - (g%south_lat+g%len_lat))* &
                         g%Rad_Earth*pi/(180.0*expdecay)))))
    enddo ; enddo
  elseif (trim(topog_config) == "halfpipe") then
    d0 = max_depth - dedge
    do i=is,ie ; do j=js,je
      d(i,j) =  dedge + d0 * abs(sin(pi*(g%geoLatT(i,j) - g%south_lat)/g%len_lat))
    enddo ; enddo
  else
    call mom_error(fatal,"initialize_topography_named: "// &
      "Unrecognized topography name "//trim(topog_config))
  endif

  ! This is here just for safety.  Hopefully it doesn't do anything.
  do i=is,ie ; do j=js,je
    if (d(i,j) > max_depth) d(i,j) = max_depth
    if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
  enddo ; enddo

  call calltree_leave(trim(mdl)//'()')
end subroutine initialize_topography_named
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
!> limit_topography ensures that  min_depth < D(x,y) < max_depth
subroutine limit_topography(D, G, param_file, max_depth, US)
  type(dyn_horgrid_type), intent(in)    :: g !< The dynamic horizontal grid type
  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
                          intent(inout) :: d !< Ocean bottom depth in m or Z if US is present
  type(param_file_type),  intent(in)    :: param_file !< Parameter file structure
  real,                   intent(in)    :: max_depth  !< Maximum depth of model in the units of D
  type(unit_scale_type), optional, intent(in) :: us   !< A dimensional unit scaling type

  ! Local variables
  real :: m_to_z  ! A dimensional rescaling factor.
  integer :: i, j
  character(len=40)  :: mdl = "limit_topography" ! This subroutine's name.
  real :: min_depth, mask_depth

  call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")

  m_to_z = 1.0 ; if (present(us)) m_to_z = us%m_to_Z

  call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
                 "If MASKING_DEPTH is unspecified, then anything shallower than "//&
                 "MINIMUM_DEPTH is assumed to be land and all fluxes are masked out. "//&
                 "If MASKING_DEPTH is specified, then all depths shallower than "//&
                 "MINIMUM_DEPTH but deeper than MASKING_DEPTH are rounded to MINIMUM_DEPTH.", &
                 units="m", default=0.0, scale=m_to_z)
  call get_param(param_file, mdl, "MASKING_DEPTH", mask_depth, &
                 "The depth below which to mask the ocean as land.", &
                 units="m", default=-9999.0, scale=m_to_z, do_not_log=.true.)

! Make sure that min_depth < D(x,y) < max_depth
  if (mask_depth < -9990.*m_to_z) then
    do j=g%jsd,g%jed ; do i=g%isd,g%ied
      d(i,j) = min( max( d(i,j), 0.5*min_depth ), max_depth )
    enddo ; enddo
  else
    do j=g%jsd,g%jed ; do i=g%isd,g%ied
      if (d(i,j)>0.) then
        d(i,j) = min( max( d(i,j), min_depth ), max_depth )
      else
        d(i,j) = 0.
      endif
    enddo ; enddo
  endif

  call calltree_leave(trim(mdl)//'()')
end subroutine limit_topography
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
!> This subroutine sets up the Coriolis parameter for a sphere
subroutine set_rotation_planetary(f, G, param_file, US)
  type(dyn_horgrid_type), intent(in)  :: g  !< The dynamic horizontal grid
  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
                          intent(out) :: f  !< Coriolis parameter (vertical component) [T-1 ~> s-1]
  type(param_file_type),  intent(in)  :: param_file !< A structure to parse for run-time parameters
  type(unit_scale_type), optional, intent(in) :: us !< A dimensional unit scaling type

! This subroutine sets up the Coriolis parameter for a sphere
  character(len=30) :: mdl = "set_rotation_planetary" ! This subroutine's name.
  integer :: i, j
  real    :: pi
  real    :: omega  ! The planetary rotation rate [T-1 ~> s-1]
  real    :: t_to_s ! A time unit conversion factor

  call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")

  t_to_s = 1.0 ; if (present(us)) t_to_s = us%T_to_s

  call get_param(param_file, "set_rotation_planetary", "OMEGA", omega, &
                 "The rotation rate of the earth.", units="s-1", &
                 default=7.2921e-5, scale=t_to_s)
  pi = 4.0*atan(1.0)

  do i=g%IsdB,g%IedB ; do j=g%JsdB,g%JedB
    f(i,j) = ( 2.0 * omega ) * sin( ( pi * g%geoLatBu(i,j) ) / 180.)
  enddo ; enddo

  call calltree_leave(trim(mdl)//'()')
end subroutine set_rotation_planetary
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
!> This subroutine sets up the Coriolis parameter for a beta-plane or f-plane
subroutine set_rotation_beta_plane(f, G, param_file, US)
  type(dyn_horgrid_type), intent(in)  :: g  !< The dynamic horizontal grid
  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
                          intent(out) :: f  !< Coriolis parameter (vertical component) [T-1 ~> s-1]
  type(param_file_type),  intent(in)  :: param_file !< A structure to parse for run-time parameters
  type(unit_scale_type), optional, intent(in) :: us !< A dimensional unit scaling type

! This subroutine sets up the Coriolis parameter for a beta-plane
  integer :: i, j
  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 m-1 ~> s-1 m-1]
  real    :: y_scl, rad_earth
  real    :: t_to_s ! A time unit conversion factor
  real    :: pi
  character(len=40)  :: mdl = "set_rotation_beta_plane" ! This subroutine's name.
  character(len=200) :: axis_units

  call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")

  t_to_s = 1.0 ; if (present(us)) t_to_s = us%T_to_s

  call get_param(param_file, mdl, "F_0", f_0, &
                 "The reference value of the Coriolis parameter with the "//&
                 "betaplane option.", units="s-1", default=0.0, scale=t_to_s)
  call get_param(param_file, mdl, "BETA", beta, &
                 "The northward gradient of the Coriolis parameter with "//&
                 "the betaplane option.", units="m-1 s-1", default=0.0, scale=t_to_s)
  call get_param(param_file, mdl, "AXIS_UNITS", axis_units, default="degrees")

  pi = 4.0*atan(1.0)
  select case (axis_units(1:1))
    case ("d")
      call get_param(param_file, mdl, "RAD_EARTH", rad_earth, &
                   "The radius of the Earth.", units="m", default=6.378e6)
      y_scl = rad_earth/pi
    case ("k"); y_scl = 1.e3
    case ("m"); y_scl = 1.
    case ("c"); y_scl = 1.e-2
    case default ; call mom_error(fatal, &
      " set_rotation_beta_plane: unknown AXIS_UNITS = "//trim(axis_units))
  end select

  do i=g%IsdB,g%IedB ; do j=g%JsdB,g%JedB
    f(i,j) = f_0 + beta * ( g%geoLatBu(i,j) * y_scl )
  enddo ; enddo

  call calltree_leave(trim(mdl)//'()')
end subroutine set_rotation_beta_plane

!> initialize_grid_rotation_angle initializes the arrays with the sine and
!!   cosine of the angle between logical north on the grid and true north.
subroutine initialize_grid_rotation_angle(G, PF)
  type(dyn_horgrid_type), intent(inout) :: g   !< The dynamic horizontal grid
  type(param_file_type),  intent(in)    :: pf  !< A structure indicating the open file
                                               !! to parse for model parameter values.

  real    :: angle, lon_scale
  real    :: len_lon    ! The periodic range of longitudes, usually 360 degrees.
  real    :: pi_720deg  ! One quarter the conversion factor from degrees to radians.
  real    :: lonb(2,2)  ! The longitude of a point, shifted to have about the same value.
  character(len=40)  :: mdl = "initialize_grid_rotation_angle" ! This subroutine's name.
  logical :: use_bugs
  integer :: i, j, m, n

  call get_param(pf, mdl, "GRID_ROTATION_ANGLE_BUGS", use_bugs, &
                 "If true, use an older algorithm to calculate the sine and "//&
                 "cosines needed rotate between grid-oriented directions and "//&
                 "true north and east.  Differences arise at the tripolar fold.", &
                 default=.false.)

  if (use_bugs) then
    do j=g%jsc,g%jec ; do i=g%isc,g%iec
      lon_scale    = cos((g%geoLatBu(i-1,j-1) + g%geoLatBu(i,j-1  ) + &
                          g%geoLatBu(i-1,j) + g%geoLatBu(i,j)) * atan(1.0)/180)
      angle        = atan2((g%geoLonBu(i-1,j) + g%geoLonBu(i,j) - &
                            g%geoLonBu(i-1,j-1) - g%geoLonBu(i,j-1))*lon_scale, &
                            g%geoLatBu(i-1,j) + g%geoLatBu(i,j) - &
                            g%geoLatBu(i-1,j-1) - g%geoLatBu(i,j-1) )
      g%sin_rot(i,j) = sin(angle) ! angle is the clockwise angle from lat/lon to ocean
      g%cos_rot(i,j) = cos(angle) ! grid (e.g. angle of ocean "north" from true north)
    enddo ; enddo

    ! This is not right at a tripolar or cubed-sphere fold.
    call pass_var(g%cos_rot, g%Domain)
    call pass_var(g%sin_rot, g%Domain)
  else
    pi_720deg = atan(1.0) / 180.0
    len_lon = 360.0 ; if (g%len_lon > 0.0) len_lon = g%len_lon
    do j=g%jsc,g%jec ; do i=g%isc,g%iec
      do n=1,2 ; do m=1,2
        lonb(m,n) = modulo_around_point(g%geoLonBu(i+m-2,j+n-2), g%geoLonT(i,j), len_lon)
      enddo ; enddo
      lon_scale = cos(pi_720deg*((g%geoLatBu(i-1,j-1) + g%geoLatBu(i,j)) + &
                                 (g%geoLatBu(i,j-1) + g%geoLatBu(i-1,j)) ) )
      angle = atan2(lon_scale*((lonb(1,2) - lonb(2,1)) + (lonb(2,2) - lonb(1,1))), &
                    (g%geoLatBu(i-1,j) - g%geoLatBu(i,j-1)) + &
                    (g%geoLatBu(i,j) - g%geoLatBu(i-1,j-1)) )
      g%sin_rot(i,j) = sin(angle) ! angle is the clockwise angle from lat/lon to ocean
      g%cos_rot(i,j) = cos(angle) ! grid (e.g. angle of ocean "north" from true north)
    enddo ; enddo

    call pass_vector(g%cos_rot, g%sin_rot, g%Domain, stagger=agrid)
  endif

end subroutine initialize_grid_rotation_angle

! -----------------------------------------------------------------------------
!> Return the modulo value of x in an interval [xc-(Lx/2) xc+(Lx/2)]
!! If Lx<=0, then it returns x without applying modulo arithmetic.
function modulo_around_point(x, xc, Lx) result(x_mod)
  real, intent(in) :: x  !< Value to which to apply modulo arithmetic
  real, intent(in) :: xc !< Center of modulo range
  real, intent(in) :: lx !< Modulo range width
  real :: x_mod          !< x shifted by an integer multiple of Lx to be close to xc.

  if (lx > 0.0) then
    x_mod = modulo(x - (xc - 0.5*lx), lx) + (xc - 0.5*lx)
  else
    x_mod = x
  endif
end function modulo_around_point

! -----------------------------------------------------------------------------
!>   This subroutine sets the open face lengths at selected points to restrict
!! passages to their observed widths based on a named set of sizes.
subroutine reset_face_lengths_named(G, param_file, name, US)
  type(dyn_horgrid_type), intent(inout) :: g  !< The dynamic horizontal grid
  type(param_file_type),  intent(in)    :: param_file !< A structure to parse for run-time parameters
  character(len=*),       intent(in)    :: name !< The name for the set of face lengths. Only "global_1deg"
                                                !! is currently implemented.
  type(unit_scale_type), optional, intent(in) :: us !< A dimensional unit scaling type

  ! Local variables
  character(len=256) :: mesg    ! Message for error messages.
  real :: m_to_l  ! A unit conversion factor [L m-1 ~> nondim]
  real :: l_to_m  ! A unit conversion factor [m L-1 ~> nondim]
  real    :: dx_2 = -1.0, dy_2 = -1.0
  real    :: pi_180
  integer :: option = -1
  integer :: i, j, isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
  pi_180 = (4.0*atan(1.0))/180.0

  select case ( trim(name) )
    case ("global_1deg")    ; option = 1 ; dx_2 = 0.5*1.0
    case default ; call mom_error(fatal, "reset_face_lengths_named: "//&
      "Unrecognized channel configuration name "//trim(name))
  end select

  m_to_l = 1.0 ; if (present(us)) m_to_l = us%m_to_L
  l_to_m = 1.0 ; if (present(us)) l_to_m = us%L_to_m

  if (option==1) then ! 1-degree settings.
    do j=jsd,jed ; do i=isdb,iedb  ! Change any u-face lengths within this loop.
      dy_2 = dx_2 * g%dyCu(i,j)*g%IdxCu(i,j) * cos(pi_180 * g%geoLatCu(i,j))

      if ((abs(g%geoLatCu(i,j)-35.5) < dy_2) .and. (g%geoLonCu(i,j) < -4.5) .and. &
          (g%geoLonCu(i,j) > -6.5)) &
        g%dy_Cu(i,j) = g%mask2dCu(i,j)*12000.0*m_to_l   ! Gibraltar

      if ((abs(g%geoLatCu(i,j)-12.5) < dy_2) .and. (abs(g%geoLonCu(i,j)-43.0) < dx_2)) &
        g%dy_Cu(i,j) = g%mask2dCu(i,j)*10000.0*m_to_l   ! Red Sea

      if ((abs(g%geoLatCu(i,j)-40.5) < dy_2) .and. (abs(g%geoLonCu(i,j)-26.0) < dx_2)) &
        g%dy_Cu(i,j) = g%mask2dCu(i,j)*5000.0*m_to_l   ! Dardanelles

      if ((abs(g%geoLatCu(i,j)-41.5) < dy_2) .and. (abs(g%geoLonCu(i,j)+220.0) < dx_2)) &
        g%dy_Cu(i,j) = g%mask2dCu(i,j)*35000.0*m_to_l   ! Tsugaru strait at 140.0e

      if ((abs(g%geoLatCu(i,j)-45.5) < dy_2) .and. (abs(g%geoLonCu(i,j)+217.5) < 0.9)) &
        g%dy_Cu(i,j) = g%mask2dCu(i,j)*15000.0*m_to_l   ! Betw Hokkaido and Sakhalin at 217&218 = 142e

      ! Greater care needs to be taken in the tripolar region.
      if ((abs(g%geoLatCu(i,j)-80.84) < 0.2) .and. (abs(g%geoLonCu(i,j)+64.9) < 0.8)) &
        g%dy_Cu(i,j) = g%mask2dCu(i,j)*38000.0*m_to_l   ! Smith Sound in Canadian Arch - tripolar region

    enddo ; enddo

    do j=jsdb,jedb ; do i=isd,ied  ! Change any v-face lengths within this loop.
      dy_2 = dx_2 * g%dyCv(i,j)*g%IdxCv(i,j) * cos(pi_180 * g%geoLatCv(i,j))
      if ((abs(g%geoLatCv(i,j)-41.0) < dy_2) .and. (abs(g%geoLonCv(i,j)-28.5) < dx_2)) &
        g%dx_Cv(i,j) = g%mask2dCv(i,j)*2500.0*m_to_l   ! Bosporus - should be 1000.0 m wide.

      if ((abs(g%geoLatCv(i,j)-13.0) < dy_2) .and. (abs(g%geoLonCv(i,j)-42.5) < dx_2)) &
        g%dx_Cv(i,j) = g%mask2dCv(i,j)*10000.0*m_to_l   ! Red Sea

      if ((abs(g%geoLatCv(i,j)+2.8) < 0.8) .and. (abs(g%geoLonCv(i,j)+241.5) < dx_2)) &
        g%dx_Cv(i,j) = g%mask2dCv(i,j)*40000.0*m_to_l   ! Makassar Straits at 241.5 W = 118.5 E

      if ((abs(g%geoLatCv(i,j)-0.56) < 0.5) .and. (abs(g%geoLonCv(i,j)+240.5) < dx_2)) &
        g%dx_Cv(i,j) = g%mask2dCv(i,j)*80000.0*m_to_l   ! entry to Makassar Straits at 240.5 W = 119.5 E

      if ((abs(g%geoLatCv(i,j)-0.19) < 0.5) .and. (abs(g%geoLonCv(i,j)+230.5) < dx_2)) &
        g%dx_Cv(i,j) = g%mask2dCv(i,j)*25000.0*m_to_l   ! Channel betw N Guinea and Halmahara 230.5 W = 129.5 E

      if ((abs(g%geoLatCv(i,j)-0.19) < 0.5) .and. (abs(g%geoLonCv(i,j)+229.5) < dx_2)) &
        g%dx_Cv(i,j) = g%mask2dCv(i,j)*25000.0*m_to_l   ! Channel betw N Guinea and Halmahara 229.5 W = 130.5 E

      if ((abs(g%geoLatCv(i,j)-0.0) < 0.25) .and. (abs(g%geoLonCv(i,j)+228.5) < dx_2)) &
        g%dx_Cv(i,j) = g%mask2dCv(i,j)*25000.0*m_to_l   ! Channel betw N Guinea and Halmahara 228.5 W = 131.5 E

      if ((abs(g%geoLatCv(i,j)+8.5) < 0.5) .and. (abs(g%geoLonCv(i,j)+244.5) < dx_2)) &
        g%dx_Cv(i,j) = g%mask2dCv(i,j)*20000.0*m_to_l   ! Lombok Straits at 244.5 W = 115.5 E

      if ((abs(g%geoLatCv(i,j)+8.5) < 0.5) .and. (abs(g%geoLonCv(i,j)+235.5) < dx_2)) &
        g%dx_Cv(i,j) = g%mask2dCv(i,j)*20000.0*m_to_l   ! Timor Straits at 235.5 W = 124.5 E

      if ((abs(g%geoLatCv(i,j)-52.5) < dy_2) .and. (abs(g%geoLonCv(i,j)+218.5) < dx_2)) &
        g%dx_Cv(i,j) = g%mask2dCv(i,j)*2500.0*m_to_l    ! Russia and Sakhalin Straits at 218.5 W = 141.5 E

      ! Greater care needs to be taken in the tripolar region.
      if ((abs(g%geoLatCv(i,j)-76.8) < 0.06) .and. (abs(g%geoLonCv(i,j)+88.7) < dx_2)) &
        g%dx_Cv(i,j) = g%mask2dCv(i,j)*8400.0*m_to_l    ! Jones Sound in Canadian Arch - tripolar region

    enddo ; enddo
  endif

  ! These checks apply regardless of the chosen option.

  do j=jsd,jed ; do i=isdb,iedb
    if (l_to_m*g%dy_Cu(i,j) > l_to_m*g%dyCu(i,j)) then
      write(mesg,'("dy_Cu of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
                   &" by ",ES11.4," at lon/lat of ", ES11.4, ES11.4)') &
                   l_to_m*g%dy_Cu(i,j), l_to_m*g%dyCu(i,j), l_to_m*g%dy_Cu(i,j)-l_to_m*g%dyCu(i,j), &
                   g%geoLonCu(i,j), g%geoLatCu(i,j)
      call mom_error(fatal,"reset_face_lengths_named "//mesg)
    endif
    g%areaCu(i,j) = g%dxCu(i,j) * g%dy_Cu(i,j)
    g%IareaCu(i,j) = 0.0
    if (g%areaCu(i,j) > 0.0) g%IareaCu(i,j) = g%mask2dCu(i,j) / (g%areaCu(i,j))
  enddo ; enddo

  do j=jsdb,jedb ; do i=isd,ied
    if (l_to_m*g%dx_Cv(i,j) > l_to_m*g%dxCv(i,j)) then
      write(mesg,'("dx_Cv of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
                   &" by ",ES11.4, " at lon/lat of ", ES11.4, ES11.4)') &
                   l_to_m*g%dx_Cv(i,j), l_to_m*g%dxCv(i,j), l_to_m*g%dx_Cv(i,j)-l_to_m*g%dxCv(i,j), &
                   g%geoLonCv(i,j), g%geoLatCv(i,j)

      call mom_error(fatal,"reset_face_lengths_named "//mesg)
    endif
    g%areaCv(i,j) = g%dyCv(i,j) * g%dx_Cv(i,j)
    g%IareaCv(i,j) = 0.0
    if (g%areaCv(i,j) > 0.0) g%IareaCv(i,j) = g%mask2dCv(i,j) / (g%areaCv(i,j))
  enddo ; enddo

end subroutine reset_face_lengths_named
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
!> This subroutine sets the open face lengths at selected points to restrict
!! passages to their observed widths from a arrays read from a file.
subroutine reset_face_lengths_file(G, param_file, US)
  type(dyn_horgrid_type), intent(inout) :: g  !< The dynamic horizontal grid
  type(param_file_type),  intent(in)    :: param_file !< A structure to parse for run-time parameters
  type(unit_scale_type), optional, intent(in) :: us !< A dimensional unit scaling type

  ! Local variables
  character(len=40)  :: mdl = "reset_face_lengths_file" ! This subroutine's name.
  character(len=256) :: mesg    ! Message for error messages.
  character(len=200) :: filename, chan_file, inputdir ! Strings for file/path
  real :: m_to_l  ! A unit conversion factor [L m-1 ~> nondim]
  real :: l_to_m  ! A unit conversion factor [m L-1 ~> nondim]
  integer :: i, j, isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
  ! These checks apply regardless of the chosen option.

  call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
  m_to_l = 1.0 ; if (present(us)) m_to_l = us%m_to_L
  l_to_m = 1.0 ; if (present(us)) l_to_m = us%L_to_m

  call get_param(param_file, mdl, "CHANNEL_WIDTH_FILE", chan_file, &
                 "The file from which the list of narrowed channels is read.", &
                 default="ocean_geometry.nc")
  call get_param(param_file,  mdl, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  filename = trim(inputdir)//trim(chan_file)
  call log_param(param_file, mdl, "INPUTDIR/CHANNEL_WIDTH_FILE", filename)

  if (is_root_pe()) then ; if (.not.file_exists(filename)) &
    call mom_error(fatal," reset_face_lengths_file: Unable to open "//&
                           trim(filename))
  endif

  call mom_read_vector(filename, "dyCuo", "dxCvo", g%dy_Cu, g%dx_Cv, g%Domain, scale=m_to_l)
  call pass_vector(g%dy_Cu, g%dx_Cv, g%Domain, to_all+scalar_pair, cgrid_ne)

  do j=jsd,jed ; do i=isdb,iedb
    if (l_to_m*g%dy_Cu(i,j) > l_to_m*g%dyCu(i,j)) then
      write(mesg,'("dy_Cu of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
                   &" by ",ES11.4," at lon/lat of ", ES11.4, ES11.4)') &
                   l_to_m*g%dy_Cu(i,j), l_to_m*g%dyCu(i,j), l_to_m*g%dy_Cu(i,j)-l_to_m*g%dyCu(i,j), &
                   g%geoLonCu(i,j), g%geoLatCu(i,j)
      call mom_error(fatal,"reset_face_lengths_file "//mesg)
    endif
    g%areaCu(i,j) = g%dxCu(i,j) * g%dy_Cu(i,j)
    g%IareaCu(i,j) = 0.0
    if (g%areaCu(i,j) > 0.0) g%IareaCu(i,j) = g%mask2dCu(i,j) / (g%areaCu(i,j))
  enddo ; enddo

  do j=jsdb,jedb ; do i=isd,ied
    if (l_to_m*g%dx_Cv(i,j) > l_to_m*g%dxCv(i,j)) then
      write(mesg,'("dx_Cv of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
                   &" by ",ES11.4, " at lon/lat of ", ES11.4, ES11.4)') &
                   l_to_m*g%dx_Cv(i,j), l_to_m*g%dxCv(i,j), l_to_m*g%dx_Cv(i,j)-l_to_m*g%dxCv(i,j), &
                   g%geoLonCv(i,j), g%geoLatCv(i,j)

      call mom_error(fatal,"reset_face_lengths_file "//mesg)
    endif
    g%areaCv(i,j) = g%dyCv(i,j) * g%dx_Cv(i,j)
    g%IareaCv(i,j) = 0.0
    if (g%areaCv(i,j) > 0.0) g%IareaCv(i,j) = g%mask2dCv(i,j) / (g%areaCv(i,j))
  enddo ; enddo

  call calltree_leave(trim(mdl)//'()')
end subroutine reset_face_lengths_file
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
!> This subroutine sets the open face lengths at selected points to restrict
!! passages to their observed widths from a list read from a file.
subroutine reset_face_lengths_list(G, param_file, US)
  type(dyn_horgrid_type), intent(inout) :: g  !< The dynamic horizontal grid
  type(param_file_type),  intent(in)    :: param_file !< A structure to parse for run-time parameters
  type(unit_scale_type), optional, intent(in)  :: us !< A dimensional unit scaling type

  ! Local variables
  character(len=120), pointer, dimension(:) :: lines => null()
  character(len=120) :: line
  character(len=200) :: filename, chan_file, inputdir, mesg ! Strings for file/path
  character(len=40)  :: mdl = "reset_face_lengths_list" ! This subroutine's name.
  real, pointer, dimension(:,:) :: &
    u_lat => null(), u_lon => null(), v_lat => null(), v_lon => null()
  real, pointer, dimension(:) :: &
    u_width => null(), v_width => null()
  real    :: m_to_l  ! A unit conversion factor [L m-1 ~> nondim]
  real    :: l_to_m  ! A unit conversion factor [m L-1 ~> nondim]
  real    :: lat, lon     ! The latitude and longitude of a point.
  real    :: len_lon      ! The periodic range of longitudes, usually 360 degrees.
  real    :: len_lat      ! The range of latitudes, usually 180 degrees.
  real    :: lon_p, lon_m ! The longitude of a point shifted by 360 degrees.
  logical :: check_360    ! If true, check for longitudes that are shifted by
                          ! +/- 360 degrees from the specified range of values.
  logical :: found_u, found_v
  logical :: unit_in_use
  integer :: ios, iounit, isu, isv
  integer :: last, num_lines, nl_read, ln, npt, u_pt, v_pt
  integer :: i, j, isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
  m_to_l = 1.0 ; if (present(us)) m_to_l = us%m_to_L
  l_to_m = 1.0 ; if (present(us)) l_to_m = us%L_to_m

  call get_param(param_file, mdl, "CHANNEL_LIST_FILE", chan_file, &
                 "The file from which the list of narrowed channels is read.", &
                 default="MOM_channel_list")
  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  filename = trim(inputdir)//trim(chan_file)
  call log_param(param_file, mdl, "INPUTDIR/CHANNEL_LIST_FILE", filename)
  call get_param(param_file, mdl, "CHANNEL_LIST_360_LON_CHECK", check_360, &
                 "If true, the channel configuration list works for any "//&
                 "longitudes in the range of -360 to 360.", default=.true.)

  if (is_root_pe()) then
    ! Open the input file.
    if (.not.file_exists(filename)) call mom_error(fatal, &
        " reset_face_lengths_list: Unable to open "//trim(filename))

    ! Find an unused unit number.
    do iounit=10,512
      INQUIRE(iounit,opened=unit_in_use) ; if (.not.unit_in_use) exit
    enddo
    if (iounit >= 512) call mom_error(fatal, &
        "reset_face_lengths_list: No unused file unit could be found.")

    ! Open the parameter file.
    open(iounit, file=trim(filename), access='SEQUENTIAL', &
         form='FORMATTED', action='READ', position='REWIND', iostat=ios)
    if (ios /= 0) call mom_error(fatal, &
            "reset_face_lengths_list: Error opening "//trim(filename))

    ! Count the number of u_width and v_width entries.
    call read_face_length_list(iounit, filename, num_lines, lines)
  endif

  len_lon = 360.0 ; if (g%len_lon > 0.0) len_lon = g%len_lon
  len_lat = 180.0 ; if (g%len_lat > 0.0) len_lat = g%len_lat
  ! Broadcast the number of lines and allocate the required space.
  call broadcast(num_lines, root_pe())
  u_pt = 0 ; v_pt = 0
  if (num_lines > 0) then
    allocate (lines(num_lines))
    if (num_lines > 0) then
      allocate(u_lat(2,num_lines)) ; u_lat(:,:) = -1e34
      allocate(u_lon(2,num_lines)) ; u_lon(:,:) = -1e34
      allocate(u_width(num_lines)) ; u_width(:) = -1e34

      allocate(v_lat(2,num_lines)) ; v_lat(:,:) = -1e34
      allocate(v_lon(2,num_lines)) ; v_lon(:,:) = -1e34
      allocate(v_width(num_lines)) ; v_width(:) = -1e34
    endif

    ! Actually read the lines.
    if (is_root_pe()) then
      call read_face_length_list(iounit, filename, nl_read, lines)
      if (nl_read /= num_lines) &
        call mom_error(fatal, 'reset_face_lengths_list : Found different '// &
                  'number of valid lines on second reading of '//trim(filename))
      close(iounit) ; iounit = -1
    endif

    ! Broadcast the lines.
    call broadcast(lines, 120, root_pe())

    ! Populate the u_width, etc., data.
    do ln=1,num_lines
      line = lines(ln)
      ! Detect keywords
      found_u = .false.; found_v = .false.
      isu = index(uppercase(line), "U_WIDTH" ); if (isu > 0) found_u = .true.
      isv = index(uppercase(line), "V_WIDTH" ); if (isv > 0) found_v = .true.

      ! Store and check the relevant values.
      if (found_u) then
        u_pt = u_pt + 1
        read(line(isu+8:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt)
        if (is_root_pe()) then
          if (check_360) then
            if ((abs(u_lon(1,u_pt)) > len_lon) .or. (abs(u_lon(2,u_pt)) > len_lon)) &
              call mom_error(warning, "reset_face_lengths_list : Out-of-bounds "//&
                 "u-longitude found when reading line "//trim(line)//" from file "//&
                 trim(filename))
            if ((abs(u_lat(1,u_pt)) > len_lat) .or. (abs(u_lat(2,u_pt)) > len_lat)) &
              call mom_error(warning, "reset_face_lengths_list : Out-of-bounds "//&
                 "u-latitude found when reading line "//trim(line)//" from file "//&
                 trim(filename))
          endif
          if (u_lat(1,u_pt) > u_lat(2,u_pt)) &
            call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
               "u-face latitudes found when reading line "//trim(line)//" from file "//&
               trim(filename))
          if (u_lon(1,u_pt) > u_lon(2,u_pt)) &
            call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
               "u-face longitudes found when reading line "//trim(line)//" from file "//&
               trim(filename))
          if (u_width(u_pt) < 0.0) &
            call mom_error(warning, "reset_face_lengths_list : Negative "//&
               "u-width found when reading line "//trim(line)//" from file "//&
               trim(filename))
        endif
      elseif (found_v) then
        v_pt = v_pt + 1
        read(line(isv+8:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt)
        if (is_root_pe()) then
          if (check_360) then
            if ((abs(v_lon(1,v_pt)) > len_lon) .or. (abs(v_lon(2,v_pt)) > len_lon)) &
              call mom_error(warning, "reset_face_lengths_list : Out-of-bounds "//&
                 "v-longitude found when reading line "//trim(line)//" from file "//&
                 trim(filename))
            if ((abs(v_lat(1,v_pt)) > len_lat) .or. (abs(v_lat(2,v_pt)) > len_lat)) &
              call mom_error(warning, "reset_face_lengths_list : Out-of-bounds "//&
                 "v-latitude found when reading line "//trim(line)//" from file "//&
                 trim(filename))
          endif
          if (v_lat(1,v_pt) > v_lat(2,v_pt)) &
            call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
               "v-face latitudes found when reading line "//trim(line)//" from file "//&
               trim(filename))
          if (v_lon(1,v_pt) > v_lon(2,v_pt)) &
            call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
               "v-face longitudes found when reading line "//trim(line)//" from file "//&
               trim(filename))
          if (v_width(v_pt) < 0.0) &
            call mom_error(warning, "reset_face_lengths_list : Negative "//&
               "v-width found when reading line "//trim(line)//" from file "//&
               trim(filename))
        endif
      endif
    enddo

    deallocate(lines)
  endif

  do j=jsd,jed ; do i=isdb,iedb
    lat = g%geoLatCu(i,j) ; lon = g%geoLonCu(i,j)
    if (check_360) then ; lon_p = lon+len_lon ; lon_m = lon-len_lon
    else ; lon_p = lon ; lon_m = lon ; endif

    do npt=1,u_pt
      if (((lat >= u_lat(1,npt)) .and. (lat <= u_lat(2,npt))) .and. &
          (((lon >= u_lon(1,npt)) .and. (lon <= u_lon(2,npt))) .or. &
           ((lon_p >= u_lon(1,npt)) .and. (lon_p <= u_lon(2,npt))) .or. &
           ((lon_m >= u_lon(1,npt)) .and. (lon_m <= u_lon(2,npt)))) ) then

        g%dy_Cu(i,j) = g%mask2dCu(i,j) * m_to_l*min(l_to_m*g%dyCu(i,j), max(u_width(npt), 0.0))
        if (j>=g%jsc .and. j<=g%jec .and. i>=g%isc .and. i<=g%iec) then ! Limit messages/checking to compute domain
          if ( g%mask2dCu(i,j) == 0.0 )  then
            write(*,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCu=0 at ",lat,lon," (",&
                u_lat(1,npt), u_lat(2,npt), u_lon(1,npt), u_lon(2,npt),") so grid metric is unmodified."
          else
            write(*,'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') &
                  "read_face_lengths_list : Modifying dy_Cu gridpoint at ",lat,lon," (",&
                  u_lat(1,npt), u_lat(2,npt), u_lon(1,npt), u_lon(2,npt),") to ",l_to_m*g%dy_Cu(i,j),"m"
          endif
        endif
      endif
    enddo

    g%areaCu(i,j) = g%dxCu(i,j) * g%dy_Cu(i,j)
    g%IareaCu(i,j) = 0.0
    if (g%areaCu(i,j) > 0.0) g%IareaCu(i,j) = g%mask2dCu(i,j) / (g%areaCu(i,j))
  enddo ; enddo

  do j=jsdb,jedb ; do i=isd,ied
    lat = g%geoLatCv(i,j) ; lon = g%geoLonCv(i,j)
    if (check_360) then ; lon_p = lon+len_lon ; lon_m = lon-len_lon
    else ; lon_p = lon ; lon_m = lon ; endif

    do npt=1,v_pt
      if (((lat >= v_lat(1,npt)) .and. (lat <= v_lat(2,npt))) .and. &
          (((lon >= v_lon(1,npt)) .and. (lon <= v_lon(2,npt))) .or. &
           ((lon_p >= v_lon(1,npt)) .and. (lon_p <= v_lon(2,npt))) .or. &
           ((lon_m >= v_lon(1,npt)) .and. (lon_m <= v_lon(2,npt)))) ) then
        g%dx_Cv(i,j) = g%mask2dCv(i,j) * m_to_l*min(l_to_m*g%dxCv(i,j), max(v_width(npt), 0.0))
        if (i>=g%isc .and. i<=g%iec .and. j>=g%jsc .and. j<=g%jec) then ! Limit messages/checking to compute domain
          if ( g%mask2dCv(i,j) == 0.0 )  then
            write(*,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCv=0 at ",lat,lon," (",&
                  v_lat(1,npt), v_lat(2,npt), v_lon(1,npt), v_lon(2,npt),") so grid metric is unmodified."
          else
            write(*,'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') &
                  "read_face_lengths_list : Modifying dx_Cv gridpoint at ",lat,lon," (",&
                  v_lat(1,npt), v_lat(2,npt), v_lon(1,npt), v_lon(2,npt),") to ",l_to_m*g%dx_Cv(i,j),"m"
          endif
        endif
      endif
    enddo

    g%areaCv(i,j) = g%dyCv(i,j) * g%dx_Cv(i,j)
    g%IareaCv(i,j) = 0.0
    if (g%areaCv(i,j) > 0.0) g%IareaCv(i,j) = g%mask2dCv(i,j) / (g%areaCv(i,j))
  enddo ; enddo

  if (num_lines > 0) then
    deallocate(u_lat) ; deallocate(u_lon) ; deallocate(u_width)
    deallocate(v_lat) ; deallocate(v_lon) ; deallocate(v_width)
  endif

  call calltree_leave(trim(mdl)//'()')
end subroutine reset_face_lengths_list
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
!>   This subroutine reads and counts the non-blank lines in the face length list file, after removing comments.
subroutine read_face_length_list(iounit, filename, num_lines, lines)
  integer,                          intent(in)  :: iounit    !< An open I/O unit number for the file
  character(len=*),                 intent(in)  :: filename  !< The name of the face-length file to read
  integer,                          intent(out) :: num_lines !< The number of non-blank lines in the file
  character(len=120), dimension(:), pointer     :: lines  !< The non-blank lines, after removing comments

  !   This subroutine reads and counts the non-blank lines in the face length
  ! list file, after removing comments.
  character(len=120) :: line, line_up
  logical :: found_u, found_v
  integer :: isu, isv, icom, verbose
  integer :: last

  num_lines = 0

  if (iounit <= 0) return
  rewind(iounit)
  do while(.true.)
    read(iounit, '(a)', end=8, err=9) line
    last = len_trim(line)
    ! Eliminate either F90 or C comments from the line.
    icom = index(line(:last), "!") ; if (icom > 0) last = icom-1
    icom = index(line(:last), "/*") ; if (icom > 0) last = icom-1
    if (last < 1) cycle

    ! Detect keywords
    line_up = uppercase(line)
    found_u = .false.; found_v = .false.
    isu = index(line_up(:last), "U_WIDTH" ); if (isu > 0) found_u = .true.
    isv = index(line_up(:last), "V_WIDTH" ); if (isv > 0) found_v = .true.

    if (found_u .and. found_v) call mom_error(fatal, &
      "read_face_length_list : both U_WIDTH and V_WIDTH found when "//&
      "reading the line "//trim(line(:last))//" in file "//trim(filename))
    if (found_u .or. found_v) then
      num_lines = num_lines + 1
      if (associated(lines)) then
        lines(num_lines) = line(1:last)
      endif
    endif
  enddo ! while (.true.)

8 continue
  return

9 call mom_error(fatal, "read_face_length_list : "//&
                  "Error while reading file "//trim(filename))

end subroutine read_face_length_list
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
!> Set the bathymetry at velocity points to be the maximum of the depths at the
!! neighoring tracer points.
subroutine set_velocity_depth_max(G)
  type(dyn_horgrid_type), intent(inout) :: g   !< The dynamic horizontal grid
  ! This subroutine sets the 4 bottom depths at velocity points to be the
  ! maximum of the adjacent depths.
  integer :: i, j

  do i=g%isd,g%ied-1 ; do j=g%jsd,g%jed
    g%Dblock_u(i,j) = g%mask2dCu(i,j) * max(g%bathyT(i,j), g%bathyT(i+1,j))
    g%Dopen_u(i,j) = g%Dblock_u(i,j)
  enddo ; enddo
  do i=g%isd,g%ied ; do j=g%jsd,g%jed-1
    g%Dblock_v(i,j) = g%mask2dCv(i,j) * max(g%bathyT(i,j), g%bathyT(i,j+1))
    g%Dopen_v(i,j) = g%Dblock_v(i,j)
  enddo ; enddo
end subroutine set_velocity_depth_max
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
!> Set the bathymetry at velocity points to be the minimum of the depths at the
!! neighoring tracer points.
subroutine set_velocity_depth_min(G)
  type(dyn_horgrid_type), intent(inout) :: g  !< The dynamic horizontal grid
  ! This subroutine sets the 4 bottom depths at velocity points to be the
  ! minimum of the adjacent depths.
  integer :: i, j

  do i=g%isd,g%ied-1 ; do j=g%jsd,g%jed
    g%Dblock_u(i,j) = g%mask2dCu(i,j) * min(g%bathyT(i,j), g%bathyT(i+1,j))
    g%Dopen_u(i,j) = g%Dblock_u(i,j)
  enddo ; enddo
  do i=g%isd,g%ied ; do j=g%jsd,g%jed-1
    g%Dblock_v(i,j) = g%mask2dCv(i,j) * min(g%bathyT(i,j), g%bathyT(i,j+1))
    g%Dopen_v(i,j) = g%Dblock_v(i,j)
  enddo ; enddo
end subroutine set_velocity_depth_min
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
!> Pre-compute global integrals of grid quantities (like masked ocean area) for
!! later use in reporting diagnostics
subroutine compute_global_grid_integrals(G, US)
  type(dyn_horgrid_type),          intent(inout) :: g  !< The dynamic horizontal grid
  type(unit_scale_type), optional, intent(in)    :: us !< A dimensional unit scaling type

  ! Local variables
  real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpforsumming
  real :: area_scale  ! A scaling factor for area into MKS units
  integer :: i,j

  area_scale = 1.0 ; if (present(us)) area_scale = us%L_to_m**2

  tmpforsumming(:,:) = 0.
  g%areaT_global = 0.0 ; g%IareaT_global = 0.0
  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    tmpforsumming(i,j) = area_scale*g%areaT(i,j) * g%mask2dT(i,j)
  enddo ; enddo
  g%areaT_global = reproducing_sum(tmpforsumming)

  if (g%areaT_global == 0.0) &
    call mom_error(fatal, "compute_global_grid_integrals: "//&
                    "zero ocean area (check topography?)")

  g%IareaT_global = 1.0 / (g%areaT_global)
end subroutine compute_global_grid_integrals
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
!> Write out a file describing the topography, Coriolis parameter, grid locations
!! and various other fixed fields from the grid.
subroutine write_ocean_geometry_file(G, param_file, directory, geom_file, US)
  type(dyn_horgrid_type),       intent(inout) :: g         !< The dynamic horizontal grid
  type(param_file_type),        intent(in)    :: param_file !< Parameter file structure
  character(len=*),             intent(in)    :: directory !< The directory into which to place the geometry file.
  character(len=*),   optional, intent(in)    :: geom_file !< If present, the name of the geometry file
                                                           !! (otherwise the file is "ocean_geometry")
  type(unit_scale_type), optional, intent(in) :: us        !< A dimensional unit scaling type

  ! Local variables.
  character(len=240) :: filepath
  character(len=40)  :: mdl = "write_ocean_geometry_file"
  integer, parameter :: nflds=23
  type(vardesc) :: vars(nflds)
  type(fieldtype) :: fields(nflds)
  real :: z_to_m_scale ! A unit conversion factor from Z to m.
  real :: s_to_t_scale ! A unit conversion factor from T-1 to s-1.
  real :: l_to_m_scale ! A unit conversion factor from L to m.
  integer :: unit
  integer :: file_threading
  integer :: nflds_used
  integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq
  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
  logical :: multiple_files
  real, dimension(G%isd :G%ied ,G%jsd :G%jed ) :: out_h
  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: out_q
  real, dimension(G%IsdB:G%IedB,G%jsd :G%jed ) :: out_u
  real, dimension(G%isd :G%ied ,G%JsdB:G%JedB) :: out_v

  call calltree_enter('write_ocean_geometry_file()')

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  z_to_m_scale = 1.0 ; if (present(us)) z_to_m_scale = us%Z_to_m
  s_to_t_scale = 1.0 ; if (present(us)) s_to_t_scale = us%s_to_T
  l_to_m_scale = 1.0 ; if (present(us)) l_to_m_scale = us%L_to_m

!   vardesc is a structure defined in MOM_io.F90.  The elements of
! this structure, in order, are:
! (1) the variable name for the NetCDF file
! (2) the variable's long name
! (3) a character indicating the  horizontal grid, which may be '1' (column),
!     'h', 'q', 'u', or 'v', for the corresponding C-grid variable
! (4) a character indicating the vertical grid, which may be 'L' (layer),
!     'i' (interface), or '1' (no vertical location)
! (5) a character indicating the time levels of the field, which may be
!    's' (snap-shot), 'p' (periodic), or '1' (no time variation)
! (6) the variable's units
  vars(1) = var_desc("geolatb","degree","latitude at corner (Bu) points",'q','1','1')
  vars(2) = var_desc("geolonb","degree","longitude at corner (Bu) points",'q','1','1')
  vars(3) = var_desc("geolat","degree", "latitude at tracer (T) points", 'h','1','1')
  vars(4) = var_desc("geolon","degree","longitude at tracer (T) points",'h','1','1')
  vars(5) = var_desc("D","meter","Basin Depth",'h','1','1')
  vars(6) = var_desc("f","s-1","Coriolis Parameter",'q','1','1')
  vars(7) = var_desc("dxCv","m","Zonal grid spacing at v points",'v','1','1')
  vars(8) = var_desc("dyCu","m","Meridional grid spacing at u points",'u','1','1')
  vars(9) = var_desc("dxCu","m","Zonal grid spacing at u points",'u','1','1')
  vars(10)= var_desc("dyCv","m","Meridional grid spacing at v points",'v','1','1')
  vars(11)= var_desc("dxT","m","Zonal grid spacing at h points",'h','1','1')
  vars(12)= var_desc("dyT","m","Meridional grid spacing at h points",'h','1','1')
  vars(13)= var_desc("dxBu","m","Zonal grid spacing at q points",'q','1','1')
  vars(14)= var_desc("dyBu","m","Meridional grid spacing at q points",'q','1','1')
  vars(15)= var_desc("Ah","m2","Area of h cells",'h','1','1')
  vars(16)= var_desc("Aq","m2","Area of q cells",'q','1','1')

  vars(17)= var_desc("dxCvo","m","Open zonal grid spacing at v points",'v','1','1')
  vars(18)= var_desc("dyCuo","m","Open meridional grid spacing at u points",'u','1','1')
  vars(19)= var_desc("wet", "nondim", "land or ocean?", 'h','1','1')

  vars(20) = var_desc("Dblock_u","m","Blocked depth at u points",'u','1','1')
  vars(21) = var_desc("Dopen_u","m","Open depth at u points",'u','1','1')
  vars(22) = var_desc("Dblock_v","m","Blocked depth at v points",'v','1','1')
  vars(23) = var_desc("Dopen_v","m","Open depth at v points",'v','1','1')


  nflds_used = 19 ; if (g%bathymetry_at_vel) nflds_used = 23

  if (present(geom_file)) then
    filepath = trim(directory) // trim(geom_file)
  else
    filepath = trim(directory) // "ocean_geometry"
  endif

  out_h(:,:) = 0.0
  out_u(:,:) = 0.0
  out_v(:,:) = 0.0
  out_q(:,:) = 0.0

  call get_param(param_file, mdl, "PARALLEL_RESTARTFILES", multiple_files, &
                 "If true, each processor writes its own restart file, "//&
                 "otherwise a single restart file is generated", &
                 default=.false.)
  file_threading = single_file
  if (multiple_files) file_threading = multiple

  call create_file(unit, trim(filepath), vars, nflds_used, fields, &
                   file_threading, dg=g)

  do j=jsq,jeq; do i=isq,ieq; out_q(i,j) = g%geoLatBu(i,j); enddo ; enddo
  call write_field(unit, fields(1), g%Domain%mpp_domain, out_q)
  do j=jsq,jeq; do i=isq,ieq; out_q(i,j) = g%geoLonBu(i,j); enddo ; enddo
  call write_field(unit, fields(2), g%Domain%mpp_domain, out_q)
  call write_field(unit, fields(3), g%Domain%mpp_domain, g%geoLatT)
  call write_field(unit, fields(4), g%Domain%mpp_domain, g%geoLonT)

  do j=js,je ; do i=is,ie ; out_h(i,j) = z_to_m_scale*g%bathyT(i,j) ; enddo ; enddo
  call write_field(unit, fields(5), g%Domain%mpp_domain, out_h)
  do j=jsq,jeq ; do i=isq,ieq ; out_q(i,j) = s_to_t_scale*g%CoriolisBu(i,j) ; enddo ; enddo
  call write_field(unit, fields(6), g%Domain%mpp_domain, out_q)

  !   I think that all of these copies are holdovers from a much earlier
  ! ancestor code in which many of the metrics were macros that could have
  ! had reduced dimensions, and that they are no longer needed in MOM6. -RWH
  do j=jsq,jeq ; do i=is,ie ; out_v(i,j) = l_to_m_scale*g%dxCv(i,j) ; enddo ; enddo
  call write_field(unit, fields(7), g%Domain%mpp_domain, out_v)
  do j=js,je ; do i=isq,ieq ; out_u(i,j) = l_to_m_scale*g%dyCu(i,j) ; enddo ; enddo
  call write_field(unit, fields(8), g%Domain%mpp_domain, out_u)

  do j=js,je ; do i=isq,ieq ; out_u(i,j) = l_to_m_scale*g%dxCu(i,j) ; enddo ; enddo
  call write_field(unit, fields(9), g%Domain%mpp_domain, out_u)
  do j=jsq,jeq ; do i=is,ie ; out_v(i,j) = l_to_m_scale*g%dyCv(i,j) ; enddo ; enddo
  call write_field(unit, fields(10), g%Domain%mpp_domain, out_v)

  do j=js,je ; do i=is,ie ; out_h(i,j) = l_to_m_scale*g%dxT(i,j); enddo ; enddo
  call write_field(unit, fields(11), g%Domain%mpp_domain, out_h)
  do j=js,je ; do i=is,ie ; out_h(i,j) = l_to_m_scale*g%dyT(i,j) ; enddo ; enddo
  call write_field(unit, fields(12), g%Domain%mpp_domain, out_h)

  do j=jsq,jeq ; do i=isq,ieq ; out_q(i,j) = l_to_m_scale*g%dxBu(i,j) ; enddo ; enddo
  call write_field(unit, fields(13), g%Domain%mpp_domain, out_q)
  do j=jsq,jeq ; do i=isq,ieq ; out_q(i,j) = l_to_m_scale*g%dyBu(i,j) ; enddo ; enddo
  call write_field(unit, fields(14), g%Domain%mpp_domain, out_q)

  do j=js,je ; do i=is,ie ; out_h(i,j) = g%areaT(i,j) ; enddo ; enddo
  call write_field(unit, fields(15), g%Domain%mpp_domain, out_h)
  do j=jsq,jeq ; do i=isq,ieq ; out_q(i,j) = g%areaBu(i,j) ; enddo ; enddo
  call write_field(unit, fields(16), g%Domain%mpp_domain, out_q)

  do j=jsq,jeq ; do i=is,ie ; out_v(i,j) = l_to_m_scale*g%dx_Cv(i,j) ; enddo ; enddo
  call write_field(unit, fields(17), g%Domain%mpp_domain, out_v)
  do j=js,je ; do i=isq,ieq ; out_u(i,j) = l_to_m_scale*g%dy_Cu(i,j) ; enddo ; enddo
  call write_field(unit, fields(18), g%Domain%mpp_domain, out_u)
  call write_field(unit, fields(19), g%Domain%mpp_domain, g%mask2dT)

  if (g%bathymetry_at_vel) then
    do j=js,je ; do i=isq,ieq ; out_u(i,j) = z_to_m_scale*g%Dblock_u(i,j) ; enddo ; enddo
    call write_field(unit, fields(20), g%Domain%mpp_domain, out_u)
    do j=js,je ; do i=isq,ieq ; out_u(i,j) = z_to_m_scale*g%Dopen_u(i,j) ; enddo ; enddo
    call write_field(unit, fields(21), g%Domain%mpp_domain, out_u)
    do j=jsq,jeq ; do i=is,ie ; out_v(i,j) = z_to_m_scale*g%Dblock_v(i,j) ; enddo ; enddo
    call write_field(unit, fields(22), g%Domain%mpp_domain, out_v)
    do j=jsq,jeq ; do i=is,ie ; out_v(i,j) = z_to_m_scale*g%Dopen_v(i,j) ; enddo ; enddo
    call write_field(unit, fields(23), g%Domain%mpp_domain, out_v)
  endif

  call close_file(unit)

  call calltree_leave('write_ocean_geometry_file()')
end subroutine write_ocean_geometry_file

end module mom_shared_initialization