MOM_grid_initialize.F90

!> Initializes horizontal grid
module mom_grid_initialize
 
! This file is part of MOM6. See LICENSE.md for the license.
 
use mom_checksums, only : hchksum, bchksum
use mom_checksums, only : uvchksum, hchksum_pair, bchksum_pair
use mom_domains, only : pass_var, pass_vector, pe_here, root_pe, broadcast
use mom_domains, only : agrid, bgrid_ne, cgrid_ne, to_all, scalar_pair
use mom_domains, only : to_north, to_south, to_east, to_west
use mom_domains, only : mom_define_domain, mom_define_io_domain
use mom_domains, only : mom_domain_type
use mom_dyn_horgrid, only : dyn_horgrid_type, set_derived_dyn_horgrid
use mom_error_handler, only : mom_error, mom_mesg, fatal, is_root_pe
use mom_error_handler, only : calltree_enter, calltree_leave
use mom_file_parser, only : get_param, log_param, log_version, param_file_type
use mom_io, only : mom_read_data, read_data, slasher, file_exists
use mom_io, only : corner, north_face, east_face
use mom_unit_scaling, only : unit_scale_type
 
use mpp_domains_mod, only : mpp_get_domain_extents, mpp_deallocate_domain
 
implicit none ; private
 
public set_grid_metrics, initialize_masks, adcroft_reciprocal
 
! 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.
 
!> Global positioning system (aka container for information to describe the grid)
type, public :: gps ; private
  real :: len_lon  !< The longitudinal or x-direction length of the domain.
  real :: len_lat  !< The latitudinal or y-direction length of the domain.
  real :: west_lon !< The western longitude of the domain or the equivalent
                   !! starting value for the x-axis.
  real :: south_lat  !< The southern latitude of the domain or the equivalent
                   !! starting value for the y-axis.
  real :: rad_earth !< The radius of the Earth [m].
  real :: lat_enhance_factor  !< The amount by which the meridional resolution
                   !! is enhanced within LAT_EQ_ENHANCE of the equator.
  real :: lat_eq_enhance !< The latitude range to the north and south of the equator
                   !! over which the resolution is enhanced, in degrees.
  logical :: isotropic !< If true, an isotropic grid on a sphere (also known as a Mercator grid)
                   !! is used. With an isotropic grid, the meridional extent of the domain
                   !! (LENLAT), the zonal extent (LENLON), and the number of grid points in each
                   !! direction are _not_ independent. In MOM the meridional extent is determined
                   !! to fit the zonal extent and the number of grid points, while grid is
                   !! perfectly isotropic.
  logical :: equator_reference !< If true, the grid is defined to have the equator at the
                   !!  nearest q- or h- grid point to (-LOWLAT*NJGLOBAL/LENLAT).
  integer :: niglobal !< The number of i-points in the global grid computational domain
  integer :: njglobal !< The number of j-points in the global grid computational domain
end type gps
 
contains
 
!> set_grid_metrics is used to set the primary values in the model's horizontal
!! grid.  The bathymetry, land-sea mask and any restricted channel widths are
!! not known yet, so these are set later.
subroutine set_grid_metrics(G, param_file, US)
  type(dyn_horgrid_type),          intent(inout) :: g  !< The dynamic horizontal grid type
  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
! This include declares and sets the variable "version".
#include "version_variable.h"
  logical :: debug
  character(len=256) :: config
 
  call calltree_enter("set_grid_metrics(), MOM_grid_initialize.F90")
  call log_version(param_file, "MOM_grid_init", version, "")
  call get_param(param_file, "MOM_grid_init", "GRID_CONFIG", config, &
                 "A character string that determines the method for "//&
                 "defining the horizontal grid.  Current options are: \n"//&
                 " \t mosaic - read the grid from a mosaic (supergrid) \n"//&
                 " \t          file set by GRID_FILE.\n"//&
                 " \t cartesian - use a (flat) Cartesian grid.\n"//&
                 " \t spherical - use a simple spherical grid.\n"//&
                 " \t mercator - use a Mercator spherical grid.", &
                 fail_if_missing=.true.)
  call get_param(param_file, "MOM_grid_init", "DEBUG", debug, &
                 "If true, write out verbose debugging data.", &
                 default=.false., debuggingparam=.true.)
 
  ! These are defaults that may be changed in the next select block.
  g%x_axis_units = "degrees_east" ; g%y_axis_units = "degrees_north"
  select case (trim(config))
    case ("mosaic");    call set_grid_metrics_from_mosaic(g, param_file, us)
    case ("cartesian"); call set_grid_metrics_cartesian(g, param_file, us)
    case ("spherical"); call set_grid_metrics_spherical(g, param_file, us)
    case ("mercator");  call set_grid_metrics_mercator(g, param_file, us)
    case ("file"); call mom_error(fatal, "MOM_grid_init: set_grid_metrics "//&
           'GRID_CONFIG "file" is no longer a supported option.  Use a '//&
           'mosaic file ("mosaic") or one of the analytic forms instead.')
    case default ; call mom_error(fatal, "MOM_grid_init: set_grid_metrics "//&
           "Unrecognized grid configuration "//trim(config))
  end select
 
  ! Calculate derived metrics (i.e. reciprocals and products)
  call calltree_enter("set_derived_metrics(), MOM_grid_initialize.F90")
  call set_derived_dyn_horgrid(g, us)
  call calltree_leave("set_derived_metrics()")
 
  if (debug) call grid_metrics_chksum('MOM_grid_init/set_grid_metrics', g, us)
 
  call calltree_leave("set_grid_metrics()")
end subroutine set_grid_metrics
 
! ------------------------------------------------------------------------------
 
!> grid_metrics_chksum performs a set of checksums on metrics on the grid for
!! debugging.
subroutine grid_metrics_chksum(parent, G, US)
  character(len=*),       intent(in) :: parent !< A string identifying the caller
  type(dyn_horgrid_type), intent(in) :: G      !< The dynamic horizontal grid type
  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
 
  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 :: halo
  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
 
  halo = min(g%ied-g%iec, g%jed-g%jec, 1)
 
  call hchksum_pair(trim(parent)//': d[xy]T', g%dxT, g%dyT, g%HI, &
                    haloshift=halo, scale=l_to_m, scalar_pair=.true.)
 
  call uvchksum(trim(parent)//': dxC[uv]', g%dxCu, g%dyCv, g%HI, haloshift=halo, scale=l_to_m)
 
  call uvchksum(trim(parent)//': dxC[uv]', g%dyCu, g%dxCv, g%HI, haloshift=halo, scale=l_to_m)
 
  call bchksum_pair(trim(parent)//': dxB[uv]', g%dxBu, g%dyBu, g%HI, haloshift=halo, scale=l_to_m)
 
  call hchksum_pair(trim(parent)//': Id[xy]T', g%IdxT, g%IdyT, g%HI, &
                    haloshift=halo, scale=m_to_l, scalar_pair=.true.)
 
  call uvchksum(trim(parent)//': Id[xy]C[uv]', g%IdxCu, g%IdyCv, g%HI, haloshift=halo, scale=m_to_l)
 
  call uvchksum(trim(parent)//': Id[xy]C[uv]', g%IdyCu, g%IdxCv, g%HI, haloshift=halo, scale=m_to_l)
 
  call bchksum_pair(trim(parent)//': Id[xy]B[uv]', g%IdxBu, g%IdyBu, g%HI, haloshift=halo, scale=m_to_l)
 
  call hchksum(g%areaT, trim(parent)//': areaT',g%HI, haloshift=halo, scale=l_to_m**2)
  call bchksum(g%areaBu, trim(parent)//': areaBu',g%HI, haloshift=halo, scale=l_to_m**2)
 
  call hchksum(g%IareaT, trim(parent)//': IareaT',g%HI, haloshift=halo, scale=m_to_l**2)
  call bchksum(g%IareaBu, trim(parent)//': IareaBu',g%HI, haloshift=halo, scale=m_to_l**2)
 
  call hchksum(g%geoLonT,trim(parent)//': geoLonT',g%HI, haloshift=halo)
  call hchksum(g%geoLatT,trim(parent)//': geoLatT',g%HI, haloshift=halo)
 
  call bchksum(g%geoLonBu, trim(parent)//': geoLonBu',g%HI, haloshift=halo)
  call bchksum(g%geoLatBu, trim(parent)//': geoLatBu',g%HI, haloshift=halo)
 
  call uvchksum(trim(parent)//': geoLonC[uv]', g%geoLonCu, g%geoLonCv, g%HI, haloshift=halo)
 
  call uvchksum(trim(parent)//': geoLatC[uv]', g%geoLatCu, g%geoLatCv, g%HI, haloshift=halo)
 
end subroutine grid_metrics_chksum
 
! ------------------------------------------------------------------------------
 
!> Sets the grid metrics from a mosaic file.
subroutine set_grid_metrics_from_mosaic(G, param_file, US)
  type(dyn_horgrid_type), intent(inout) :: G           !< The dynamic horizontal grid type
  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, dimension(G%isd :G%ied ,G%jsd :G%jed ) :: tempH1, tempH2, tempH3, tempH4
  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: tempQ1, tempQ2, tempQ3, tempQ4
  real, dimension(G%IsdB:G%IedB,G%jsd :G%jed ) :: tempE1, tempE2
  real, dimension(G%isd :G%ied ,G%JsdB:G%JedB) :: tempN1, tempN2
  ! These arrays are a holdover from earlier code in which the arrays in G were
  ! macros and may have had reduced dimensions.
  real, dimension(G%isd :G%ied ,G%jsd :G%jed ) :: dxT, dyT, areaT
  real, dimension(G%IsdB:G%IedB,G%jsd :G%jed ) :: dxCu, dyCu
  real, dimension(G%isd :G%ied ,G%JsdB:G%JedB) :: dxCv, dyCv
  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: dxBu, dyBu, areaBu
  ! This are symmetric arrays, corresponding to the data in the mosaic file
  real, dimension(2*G%isd-2:2*G%ied+1,2*G%jsd-2:2*G%jed+1) :: tmpT
  real, dimension(2*G%isd-3:2*G%ied+1,2*G%jsd-2:2*G%jed+1) :: tmpU
  real, dimension(2*G%isd-2:2*G%ied+1,2*G%jsd-3:2*G%jed+1) :: tmpV
  real, dimension(2*G%isd-3:2*G%ied+1,2*G%jsd-3:2*G%jed+1) :: tmpZ
  real, dimension(:,:), allocatable :: tmpGlbl
  real :: m_to_L  ! A unit conversion factor [L m-1 ~> nondim]
  character(len=200) :: filename, grid_file, inputdir
  character(len=64)  :: mdl = "MOM_grid_init set_grid_metrics_from_mosaic"
  integer :: err=0, ni, nj, global_indices(4)
  type(mom_domain_type) :: SGdom ! Supergrid domain
  logical :: lon_bug  ! If true use an older buggy answer in the tripolar longitude.
  integer :: i, j, i2, j2
  integer :: npei,npej
  integer, dimension(:), allocatable :: exni,exnj
  integer        :: start(4), nread(4)
 
  call calltree_enter("set_grid_metrics_from_mosaic(), MOM_grid_initialize.F90")
 
  m_to_l = 1.0 ; if (present(us)) m_to_l = us%m_to_L
  call get_param(param_file, mdl, "GRID_FILE", grid_file, &
                 "Name of the file from which to read horizontal grid data.", &
                 fail_if_missing=.true.)
  call get_param(param_file, mdl, "USE_TRIPOLAR_GEOLONB_BUG", lon_bug, &
                 "If true, use older code that incorrectly sets the longitude "//&
                 "in some points along the tripolar fold to be off by 360 degrees.", &
                 default=.false.)
  call get_param(param_file,  mdl, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  filename = trim(adjustl(inputdir)) // trim(adjustl(grid_file))
  call log_param(param_file, mdl, "INPUTDIR/GRID_FILE", filename)
  if (.not.file_exists(filename)) &
    call mom_error(fatal," set_grid_metrics_from_mosaic: Unable to open "//&
                           trim(filename))
 
  ! Initialize everything to 0.
  dxcu(:,:) = 0.0 ; dycu(:,:) = 0.0
  dxcv(:,:) = 0.0 ; dycv(:,:) = 0.0
  dxbu(:,:) = 0.0 ; dybu(:,:) = 0.0 ; areabu(:,:) = 0.0
 
  !<MISSING CODE TO READ REFINEMENT LEVEL>
  ni = 2*(g%iec-g%isc+1) ! i size of supergrid
  nj = 2*(g%jec-g%jsc+1) ! j size of supergrid
 
  ! Define a domain for the supergrid (SGdom)
  npei = g%domain%layout(1) ; npej = g%domain%layout(2)
  allocate(exni(npei)) ; allocate(exnj(npej))
  call mpp_get_domain_extents(g%domain%mpp_domain, exni, exnj)
  allocate(sgdom%mpp_domain)
  sgdom%nihalo = 2*g%domain%nihalo+1
  sgdom%njhalo = 2*g%domain%njhalo+1
  sgdom%niglobal = 2*g%domain%niglobal
  sgdom%njglobal = 2*g%domain%njglobal
  sgdom%layout(:) = g%domain%layout(:)
  sgdom%io_layout(:) = g%domain%io_layout(:)
  global_indices(1) = 1+sgdom%nihalo
  global_indices(2) = sgdom%niglobal+sgdom%nihalo
  global_indices(3) = 1+sgdom%njhalo
  global_indices(4) = sgdom%njglobal+sgdom%njhalo
  exni(:) = 2*exni(:) ; exnj(:) = 2*exnj(:)
  if (associated(g%domain%maskmap)) then
     call mom_define_domain(global_indices, sgdom%layout, sgdom%mpp_domain, &
            xflags=g%domain%X_FLAGS, yflags=g%domain%Y_FLAGS, &
            xhalo=sgdom%nihalo, yhalo=sgdom%njhalo, &
            xextent=exni,yextent=exnj, &
            symmetry=.true., name="MOM_MOSAIC", maskmap=g%domain%maskmap)
  else
     call mom_define_domain(global_indices, sgdom%layout, sgdom%mpp_domain, &
            xflags=g%domain%X_FLAGS, yflags=g%domain%Y_FLAGS, &
            xhalo=sgdom%nihalo, yhalo=sgdom%njhalo, &
            xextent=exni,yextent=exnj, &
            symmetry=.true., name="MOM_MOSAIC")
  endif
 
  call mom_define_io_domain(sgdom%mpp_domain, sgdom%io_layout)
  deallocate(exni)
  deallocate(exnj)
 
  ! Read X from the supergrid
  tmpz(:,:) = 999.
  call mom_read_data(filename, 'x', tmpz, sgdom, position=corner)
 
  if (lon_bug) then
    call pass_var(tmpz, sgdom, position=corner)
  else
    call pass_var(tmpz, sgdom, position=corner, inner_halo=0)
  endif
  call extrapolate_metric(tmpz, 2*(g%jsc-g%jsd)+2, missing=999.)
  do j=g%jsd,g%jed ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
    g%geoLonT(i,j) = tmpz(i2-1,j2-1)
  enddo ; enddo
  do j=g%JsdB,g%JedB ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
    g%geoLonBu(i,j) = tmpz(i2,j2)
  enddo ; enddo
  do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
    g%geoLonCu(i,j) = tmpz(i2,j2-1)
  enddo ; enddo
  do j=g%JsdB,g%JedB ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
    g%geoLonCv(i,j) = tmpz(i2-1,j2)
  enddo ; enddo
  ! For some reason, this messes up the solution...
  !   call pass_var(G%geoLonBu, G%domain, position=CORNER)
 
  ! Read Y from the supergrid
  tmpz(:,:) = 999.
  call mom_read_data(filename, 'y', tmpz, sgdom, position=corner)
 
  call pass_var(tmpz, sgdom, position=corner)
  call extrapolate_metric(tmpz, 2*(g%jsc-g%jsd)+2, missing=999.)
  do j=g%jsd,g%jed ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
    g%geoLatT(i,j) = tmpz(i2-1,j2-1)
  enddo ; enddo
  do j=g%JsdB,g%JedB ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
    g%geoLatBu(i,j) = tmpz(i2,j2)
  enddo ; enddo
  do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
    g%geoLatCu(i,j) = tmpz(i2,j2-1)
  enddo ; enddo
  do j=g%JsdB,g%JedB ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
    g%geoLatCv(i,j) = tmpz(i2-1,j2)
  enddo ; enddo
 
  ! Read DX,DY from the supergrid
  tmpu(:,:) = 0. ; tmpv(:,:) = 0.
  call mom_read_data(filename,'dx',tmpv,sgdom,position=north_face)
  call mom_read_data(filename,'dy',tmpu,sgdom,position=east_face)
  call pass_vector(tmpu, tmpv, sgdom, to_all+scalar_pair, cgrid_ne)
  call extrapolate_metric(tmpv, 2*(g%jsc-g%jsd)+2, missing=0.)
  call extrapolate_metric(tmpu, 2*(g%jsc-g%jsd)+2, missing=0.)
 
  do j=g%jsd,g%jed ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
    dxt(i,j) = tmpv(i2-1,j2-1) + tmpv(i2,j2-1)
    dyt(i,j) = tmpu(i2-1,j2-1) + tmpu(i2-1,j2)
  enddo ; enddo
 
  do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
    dxcu(i,j) = tmpv(i2,j2-1) + tmpv(i2+1,j2-1)
    dycu(i,j) = tmpu(i2,j2-1) + tmpu(i2,j2)
  enddo ; enddo
 
  do j=g%JsdB,g%JedB ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
    dxcv(i,j) = tmpv(i2-1,j2) + tmpv(i2,j2)
    dycv(i,j) = tmpu(i2-1,j2) + tmpu(i2-1,j2+1)
  enddo ; enddo
 
  do j=g%JsdB,g%JedB ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
    dxbu(i,j) = tmpv(i2,j2) + tmpv(i2+1,j2)
    dybu(i,j) = tmpu(i2,j2) + tmpu(i2,j2+1)
  enddo ; enddo
 
  ! Read AREA from the supergrid
  tmpt(:,:) = 0.
  call mom_read_data(filename, 'area', tmpt, sgdom)
  call pass_var(tmpt, sgdom)
  call extrapolate_metric(tmpt, 2*(g%jsc-g%jsd)+2, missing=0.)
 
  do j=g%jsd,g%jed ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
    areat(i,j) = (tmpt(i2-1,j2-1) + tmpt(i2,j2)) + &
                 (tmpt(i2-1,j2) + tmpt(i2,j2-1))
  enddo ; enddo
  do j=g%JsdB,g%JedB ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
    areabu(i,j) = (tmpt(i2,j2) + tmpt(i2+1,j2+1)) + &
                  (tmpt(i2,j2+1) + tmpt(i2+1,j2))
  enddo ; enddo
 
  ni=sgdom%niglobal
  nj=sgdom%njglobal
  call mpp_deallocate_domain(sgdom%mpp_domain)
  deallocate(sgdom%mpp_domain)
 
  call pass_vector(dycu, dxcv, g%Domain, to_all+scalar_pair, cgrid_ne)
  call pass_vector(dxcu, dycv, g%Domain, to_all+scalar_pair, cgrid_ne)
  call pass_vector(dxbu, dybu, g%Domain, to_all+scalar_pair, bgrid_ne)
  call pass_var(areat, g%Domain)
  call pass_var(areabu, g%Domain, position=corner)
 
  do i=g%isd,g%ied ; do j=g%jsd,g%jed
    g%dxT(i,j) = m_to_l*dxt(i,j) ; g%dyT(i,j) = m_to_l*dyt(i,j) ; g%areaT(i,j) = m_to_l**2*areat(i,j)
  enddo ; enddo
  do i=g%IsdB,g%IedB ; do j=g%jsd,g%jed
    g%dxCu(i,j) = m_to_l*dxcu(i,j) ; g%dyCu(i,j) = m_to_l*dycu(i,j)
  enddo ; enddo
  do i=g%isd,g%ied ; do j=g%JsdB,g%JedB
    g%dxCv(i,j) = m_to_l*dxcv(i,j) ; g%dyCv(i,j) = m_to_l*dycv(i,j)
  enddo ; enddo
  do i=g%IsdB,g%IedB ; do j=g%JsdB,g%JedB
    g%dxBu(i,j) = m_to_l*dxbu(i,j) ; g%dyBu(i,j) = m_to_l*dybu(i,j) ; g%areaBu(i,j) = m_to_l**2*areabu(i,j)
  enddo ; enddo
 
  ! Construct axes for diagnostic output (only necessary because "ferret" uses
  ! broken convention for interpretting netCDF files).
  start(:) = 1 ; nread(:) = 1
  start(2) = 2 ; nread(1) = ni+1 ; nread(2) = 2
  allocate( tmpglbl(ni+1,2) )
  if (is_root_pe()) &
    call read_data(filename, "x", tmpglbl, start, nread, no_domain=.true.)
  call broadcast(tmpglbl, 2*(ni+1), root_pe())
 
  ! I don't know why the second axis is 1 or 2 here. -RWH
  do i=g%isg,g%ieg
    g%gridLonT(i) = tmpglbl(2*(i-g%isg)+2,2)
  enddo
  ! Note that the dynamic grid always uses symmetric memory for the global
  ! arrays G%gridLatB and G%gridLonB.
  do i=g%isg-1,g%ieg
    g%gridLonB(i) = tmpglbl(2*(i-g%isg)+3,1)
  enddo
  deallocate( tmpglbl )
 
  allocate( tmpglbl(1, nj+1) )
  start(:) = 1 ; nread(:) = 1
  start(1) = int(ni/4)+1 ; nread(2) = nj+1
  if (is_root_pe()) &
    call read_data(filename, "y", tmpglbl, start, nread, no_domain=.true.)
  call broadcast(tmpglbl, nj+1, root_pe())
 
  do j=g%jsg,g%jeg
    g%gridLatT(j) = tmpglbl(1,2*(j-g%jsg)+2)
  enddo
  do j=g%jsg-1,g%jeg
    g%gridLatB(j) = tmpglbl(1,2*(j-g%jsg)+3)
  enddo
  deallocate( tmpglbl )
 
  call calltree_leave("set_grid_metrics_from_mosaic()")
end subroutine set_grid_metrics_from_mosaic
 
 
! ------------------------------------------------------------------------------
 
!> Calculate the values of the metric terms for a Cartesian grid that
!! might be used and save them in arrays.
!!
!! Within this subroutine, the x- and y- grid spacings and their
!! inverses and the cell areas centered on h, q, u, and v points are
!! calculated, as are the geographic locations of each of these 4
!! sets of points.
subroutine set_grid_metrics_cartesian(G, param_file, US)
  type(dyn_horgrid_type), intent(inout) :: G           !< The dynamic horizontal grid type
  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
  integer :: i, j, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, I1off, J1off
  integer :: niglobal, njglobal
  real :: grid_latT(G%jsd:G%jed), grid_latB(G%JsdB:G%JedB)
  real :: grid_lonT(G%isd:G%ied), grid_lonB(G%IsdB:G%IedB)
  real :: dx_everywhere, dy_everywhere ! Grid spacings [m].
  real :: I_dx, I_dy                   ! Inverse grid spacings [m-1].
  real :: PI
  real :: m_to_L  ! A unit conversion factor [L m-1 ~> nondim]
  real :: L_to_m  ! A unit conversion factor [m L-1 ~> nondim]
  character(len=80) :: units_temp
  character(len=48) :: mdl  = "MOM_grid_init set_grid_metrics_cartesian"
 
  niglobal = g%Domain%niglobal ; njglobal = g%Domain%njglobal
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
  i1off = g%idg_offset ; j1off = g%jdg_offset
 
  call calltree_enter("set_grid_metrics_cartesian(), MOM_grid_initialize.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
  pi = 4.0*atan(1.0)
 
  call get_param(param_file, mdl, "AXIS_UNITS", units_temp, &
                 "The units for the Cartesian axes. Valid entries are: \n"//&
                 " \t degrees - degrees of latitude and longitude \n"//&
                 " \t m - meters \n \t k - kilometers", default="degrees")
  call get_param(param_file, mdl, "SOUTHLAT", g%south_lat, &
                 "The southern latitude of the domain or the equivalent "//&
                 "starting value for the y-axis.", units=units_temp, &
                 fail_if_missing=.true.)
  call get_param(param_file, mdl, "LENLAT", g%len_lat, &
                 "The latitudinal or y-direction length of the domain.", &
                 units=units_temp, fail_if_missing=.true.)
  call get_param(param_file, mdl, "WESTLON", g%west_lon, &
                 "The western longitude of the domain or the equivalent "//&
                 "starting value for the x-axis.", units=units_temp, &
                 default=0.0)
  call get_param(param_file, mdl, "LENLON", g%len_lon, &
                 "The longitudinal or x-direction length of the domain.", &
                 units=units_temp, fail_if_missing=.true.)
  call get_param(param_file, mdl, "RAD_EARTH", g%Rad_Earth, &
                 "The radius of the Earth.", units="m", default=6.378e6)
 
  if (units_temp(1:1) == 'k') then
    g%x_axis_units = "kilometers" ; g%y_axis_units = "kilometers"
  elseif (units_temp(1:1) == 'm') then
    g%x_axis_units = "meters" ; g%y_axis_units = "meters"
  endif
  call log_param(param_file, mdl, "explicit AXIS_UNITS", g%x_axis_units)
 
  ! Note that the dynamic grid always uses symmetric memory for the global
  ! arrays G%gridLatB and G%gridLonB.
  do j=g%jsg-1,g%jeg
    g%gridLatB(j) = g%south_lat + g%len_lat*real(j-(g%jsg-1))/real(njglobal)
  enddo
  do j=g%jsg,g%jeg
    g%gridLatT(j) = g%south_lat + g%len_lat*(real(j-g%jsg)+0.5)/real(njglobal)
  enddo
  do i=g%isg-1,g%ieg
    g%gridLonB(i) = g%west_lon + g%len_lon*real(i-(g%isg-1))/real(niglobal)
  enddo
  do i=g%isg,g%ieg
    g%gridLonT(i) = g%west_lon + g%len_lon*(real(i-g%isg)+0.5)/real(niglobal)
  enddo
 
  do j=jsdb,jedb
    grid_latb(j) = g%south_lat + g%len_lat*real(j+j1off-(g%jsg-1))/real(njglobal)
  enddo
  do j=jsd,jed
    grid_latt(j) = g%south_lat + g%len_lat*(real(j+j1off-g%jsg)+0.5)/real(njglobal)
  enddo
  do i=isdb,iedb
    grid_lonb(i) = g%west_lon + g%len_lon*real(i+i1off-(g%isg-1))/real(niglobal)
  enddo
  do i=isd,ied
    grid_lont(i) = g%west_lon + g%len_lon*(real(i+i1off-g%isg)+0.5)/real(niglobal)
  enddo
 
  if (units_temp(1:1) == 'k') then ! Axes are measured in km.
    dx_everywhere = 1000.0 * g%len_lon / (real(niglobal))
    dy_everywhere = 1000.0 * g%len_lat / (real(njglobal))
  elseif (units_temp(1:1) == 'm') then ! Axes are measured in m.
    dx_everywhere = g%len_lon / (real(niglobal))
    dy_everywhere = g%len_lat / (real(njglobal))
  else ! Axes are measured in degrees of latitude and longitude.
    dx_everywhere = g%Rad_Earth * g%len_lon * pi / (180.0 * niglobal)
    dy_everywhere = g%Rad_Earth * g%len_lat * pi / (180.0 * njglobal)
  endif
 
  i_dx = 1.0 / dx_everywhere ; i_dy = 1.0 / dy_everywhere
 
  do j=jsdb,jedb ; do i=isdb,iedb
    g%geoLonBu(i,j) = grid_lonb(i) ; g%geoLatBu(i,j) = grid_latb(j)
 
    g%dxBu(i,j) = m_to_l*dx_everywhere ; g%IdxBu(i,j) = l_to_m*i_dx
    g%dyBu(i,j) = m_to_l*dy_everywhere ; g%IdyBu(i,j) = l_to_m*i_dy
    g%areaBu(i,j) = m_to_l**2*dx_everywhere * dy_everywhere ; g%IareaBu(i,j) = l_to_m**2*i_dx * i_dy
  enddo ; enddo
 
  do j=jsd,jed ; do i=isd,ied
    g%geoLonT(i,j) = grid_lont(i) ; g%geoLatT(i,j) = grid_latt(j)
    g%dxT(i,j) = m_to_l*dx_everywhere ; g%IdxT(i,j) = l_to_m*i_dx
    g%dyT(i,j) = m_to_l*dy_everywhere ; g%IdyT(i,j) = l_to_m*i_dy
    g%areaT(i,j) = m_to_l**2*dx_everywhere * dy_everywhere ; g%IareaT(i,j) = l_to_m**2*i_dx * i_dy
  enddo ; enddo
 
  do j=jsd,jed ; do i=isdb,iedb
    g%geoLonCu(i,j) = grid_lonb(i) ; g%geoLatCu(i,j) = grid_latt(j)
 
    g%dxCu(i,j) = m_to_l*dx_everywhere ; g%IdxCu(i,j) = l_to_m*i_dx
    g%dyCu(i,j) = m_to_l*dy_everywhere ; g%IdyCu(i,j) = l_to_m*i_dy
  enddo ; enddo
 
  do j=jsdb,jedb ; do i=isd,ied
    g%geoLonCv(i,j) = grid_lont(i) ; g%geoLatCv(i,j) = grid_latb(j)
 
    g%dxCv(i,j) = m_to_l*dx_everywhere ; g%IdxCv(i,j) = l_to_m*i_dx
    g%dyCv(i,j) = m_to_l*dy_everywhere ; g%IdyCv(i,j) = l_to_m*i_dy
  enddo ; enddo
 
  call calltree_leave("set_grid_metrics_cartesian()")
end subroutine set_grid_metrics_cartesian
 
! ------------------------------------------------------------------------------
 
!> Calculate the values of the metric terms that might be used
!! and save them in arrays.
!!
!! Within this subroutine, the x- and y- grid spacings and their
!! inverses and the cell areas centered on h, q, u, and v points are
!! calculated, as are the geographic locations of each of these 4
!! sets of points.
subroutine set_grid_metrics_spherical(G, param_file, US)
  type(dyn_horgrid_type), intent(inout) :: G           !< The dynamic horizontal grid type
  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 :: PI, PI_180! PI = 3.1415926... as 4*atan(1)
  integer :: i, j, isd, ied, jsd, jed
  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, IsdB, IedB, JsdB, JedB
  integer :: i_offset, j_offset
  real :: grid_latT(G%jsd:G%jed), grid_latB(G%JsdB:G%JedB)
  real :: grid_lonT(G%isd:G%ied), grid_lonB(G%IsdB:G%IedB)
  real :: dLon,dLat,latitude,longitude,dL_di
  real :: m_to_L  ! A unit conversion factor [L m-1 ~> nondim]
  character(len=48)  :: mdl  = "MOM_grid_init set_grid_metrics_spherical"
 
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
  i_offset = g%idg_offset ; j_offset = g%jdg_offset
 
  call calltree_enter("set_grid_metrics_spherical(), MOM_grid_initialize.F90")
  m_to_l = 1.0 ; if (present(us)) m_to_l = us%m_to_L
 
!    Calculate the values of the metric terms that might be used
!  and save them in arrays.
  pi = 4.0*atan(1.0); pi_180 = atan(1.0)/45.
 
  call get_param(param_file, mdl, "SOUTHLAT", g%south_lat, &
                 "The southern latitude of the domain.", units="degrees", &
                 fail_if_missing=.true.)
  call get_param(param_file, mdl, "LENLAT", g%len_lat, &
                 "The latitudinal length of the domain.", units="degrees", &
                 fail_if_missing=.true.)
  call get_param(param_file, mdl, "WESTLON", g%west_lon, &
                 "The western longitude of the domain.", units="degrees", &
                 default=0.0)
  call get_param(param_file, mdl, "LENLON", g%len_lon, &
                 "The longitudinal length of the domain.", units="degrees", &
                 fail_if_missing=.true.)
  call get_param(param_file, mdl, "RAD_EARTH", g%Rad_Earth, &
                 "The radius of the Earth.", units="m", default=6.378e6)
 
  dlon = g%len_lon/g%Domain%niglobal
  dlat = g%len_lat/g%Domain%njglobal
 
  ! Note that the dynamic grid always uses symmetric memory for the global
  ! arrays G%gridLatB and G%gridLonB.
  do j=g%jsg-1,g%jeg
    latitude = g%south_lat + dlat*(real(j-(g%jsg-1)))
    g%gridLatB(j) = min(max(latitude,-90.),90.)
  enddo
  do j=g%jsg,g%jeg
    latitude = g%south_lat + dlat*(real(j-g%jsg)+0.5)
    g%gridLatT(j) = min(max(latitude,-90.),90.)
  enddo
  do i=g%isg-1,g%ieg
    g%gridLonB(i) = g%west_lon + dlon*(real(i-(g%isg-1)))
  enddo
  do i=g%isg,g%ieg
    g%gridLonT(i) = g%west_lon + dlon*(real(i-g%isg)+0.5)
  enddo
 
  do j=jsdb,jedb
    latitude = g%south_lat + dlat* real(j+j_offset-(g%jsg-1))
    grid_latb(j) = min(max(latitude,-90.),90.)
  enddo
  do j=jsd,jed
    latitude = g%south_lat + dlat*(real(j+j_offset-g%jsg)+0.5)
    grid_latt(j) = min(max(latitude,-90.),90.)
  enddo
  do i=isdb,iedb
    grid_lonb(i) = g%west_lon + dlon*real(i+i_offset-(g%isg-1))
  enddo
  do i=isd,ied
    grid_lont(i) = g%west_lon + dlon*(real(i+i_offset-g%isg)+0.5)
  enddo
 
  dl_di = (g%len_lon * 4.0*atan(1.0)) / (180.0 * g%Domain%niglobal)
  do j=jsdb,jedb ; do i=isdb,iedb
    g%geoLonBu(i,j) = grid_lonb(i)
    g%geoLatBu(i,j) = grid_latb(j)
 
    ! The following line is needed to reproduce the solution from
    ! set_grid_metrics_mercator when used to generate a simple spherical grid.
    g%dxBu(i,j) = m_to_l*g%Rad_Earth * cos( g%geoLatBu(i,j)*pi_180 ) * dl_di
!   G%dxBu(I,J) = m_to_L*G%Rad_Earth * dLon*PI_180 * COS( G%geoLatBu(I,J)*PI_180 )
    g%dyBu(i,j) = m_to_l*g%Rad_Earth * dlat*pi_180
    g%areaBu(i,j) = g%dxBu(i,j) * g%dyBu(i,j)
  enddo ; enddo
 
  do j=jsdb,jedb ; do i=isd,ied
    g%geoLonCv(i,j) = grid_lont(i)
    g%geoLatCv(i,j) = grid_latb(j)
 
    ! The following line is needed to reproduce the solution from
    ! set_grid_metrics_mercator when used to generate a simple spherical grid.
    g%dxCv(i,j) = m_to_l*g%Rad_Earth * cos( g%geoLatCv(i,j)*pi_180 ) * dl_di
!   G%dxCv(i,J) = m_to_L*G%Rad_Earth * (dLon*PI_180) * COS( G%geoLatCv(i,J)*PI_180 )
    g%dyCv(i,j) = m_to_l*g%Rad_Earth * dlat*pi_180
  enddo ; enddo
 
  do j=jsd,jed ; do i=isdb,iedb
    g%geoLonCu(i,j) = grid_lonb(i)
    g%geoLatCu(i,j) = grid_latt(j)
 
    ! The following line is needed to reproduce the solution from
    ! set_grid_metrics_mercator when used to generate a simple spherical grid.
    g%dxCu(i,j) = m_to_l*g%Rad_Earth * cos( g%geoLatCu(i,j)*pi_180 ) * dl_di
!   G%dxCu(I,j) = m_to_L*G%Rad_Earth * dLon*PI_180 * COS( latitude )
    g%dyCu(i,j) = m_to_l*g%Rad_Earth * dlat*pi_180
  enddo ; enddo
 
  do j=jsd,jed ; do i=isd,ied
    g%geoLonT(i,j) = grid_lont(i)
    g%geoLatT(i,j) = grid_latt(j)
 
    ! The following line is needed to reproduce the solution from
    ! set_grid_metrics_mercator when used to generate a simple spherical grid.
    g%dxT(i,j) = m_to_l*g%Rad_Earth * cos( g%geoLatT(i,j)*pi_180 ) * dl_di
!   G%dxT(i,j) = G%Rad_Earth * dLon*PI_180 * COS( latitude )
    g%dyT(i,j) = m_to_l*g%Rad_Earth * dlat*pi_180
 
!   latitude = G%geoLatCv(i,J)*PI_180             ! In radians
!   dL_di    = G%geoLatCv(i,max(jsd,J-1))*PI_180  ! In radians
!   G%areaT(i,j) = m_to_L**2 * Rad_Earth**2*dLon*dLat*ABS(SIN(latitude)-SIN(dL_di))
    g%areaT(i,j) = g%dxT(i,j) * g%dyT(i,j)
  enddo ; enddo
 
  call calltree_leave("set_grid_metrics_spherical()")
end subroutine set_grid_metrics_spherical
 
!> Calculate the values of the metric terms that might be used
!! and save them in arrays.
!!
!! Within this subroutine, the x- and y- grid spacings and their
!! inverses and the cell areas centered on h, q, u, and v points are
!! calculated, as are the geographic locations of each of these 4
!! sets of points.
subroutine set_grid_metrics_mercator(G, param_file, US)
  type(dyn_horgrid_type), intent(inout) :: G           !< The dynamic horizontal grid type
  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
  integer :: i, j, isd, ied, jsd, jed
  integer :: I_off, J_off
  type(gps) :: GP
  character(len=128) :: warnmesg
  character(len=48)  :: mdl = "MOM_grid_init set_grid_metrics_mercator"
  real :: PI, PI_2! PI = 3.1415926... as 4*atan(1), PI_2 = (PI) /2.0
  real :: y_q, y_h, jd, x_q, x_h, id
  real, dimension(G%isd:G%ied,G%jsd:G%jed) :: &
    xh, yh ! Latitude and longitude of h points in radians.
  real, dimension(G%IsdB:G%IedB,G%jsd:G%jed) :: &
    xu, yu ! Latitude and longitude of u points in radians.
  real, dimension(G%isd:G%ied,G%JsdB:G%JedB) :: &
    xv, yv ! Latitude and longitude of v points in radians.
  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: &
    xq, yq ! Latitude and longitude of q points in radians.
  real :: fnRef           ! fnRef is the value of Int_dj_dy or
                          ! Int_dj_dy at a latitude or longitude that is
  real :: jRef, iRef      ! being set to be at grid index jRef or iRef.
  real :: m_to_L  ! A unit conversion factor [L m-1 ~> nondim]
  integer :: itt1, itt2
  logical :: debug = .false., simple_area = .true.
  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, IsdB, IedB, JsdB, JedB
 
  !   All of the metric terms should be defined over the domain from
  ! isd to ied.  Outside of the physical domain, both the metrics
  ! and their inverses may be set to zero.
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
  i_off = g%idg_offset ; j_off = g%jdg_offset
 
  gp%niglobal = g%Domain%niglobal
  gp%njglobal = g%Domain%njglobal
 
  call calltree_enter("set_grid_metrics_mercator(), MOM_grid_initialize.F90")
 
  m_to_l = 1.0 ; if (present(us)) m_to_l = us%m_to_L
  !   Calculate the values of the metric terms that might be used
  ! and save them in arrays.
  pi = 4.0*atan(1.0) ; pi_2 = 0.5*pi
 
  call get_param(param_file, mdl, "SOUTHLAT", gp%south_lat, &
                 "The southern latitude of the domain.", units="degrees", &
                 fail_if_missing=.true.)
  call get_param(param_file, mdl, "LENLAT", gp%len_lat, &
                 "The latitudinal length of the domain.", units="degrees", &
                 fail_if_missing=.true.)
  call get_param(param_file, mdl, "WESTLON", gp%west_lon, &
                 "The western longitude of the domain.", units="degrees", &
                 default=0.0)
  call get_param(param_file, mdl, "LENLON", gp%len_lon, &
                 "The longitudinal length of the domain.", units="degrees", &
                 fail_if_missing=.true.)
  call get_param(param_file, mdl, "RAD_EARTH", gp%Rad_Earth, &
                 "The radius of the Earth.", units="m", default=6.378e6)
  g%south_lat = gp%south_lat ; g%len_lat = gp%len_lat
  g%west_lon = gp%west_lon ; g%len_lon = gp%len_lon
  g%Rad_Earth = gp%Rad_Earth
  call get_param(param_file, mdl, "ISOTROPIC", gp%isotropic, &
                 "If true, an isotropic grid on a sphere (also known as "//&
                 "a Mercator grid) is used. With an isotropic grid, the "//&
                 "meridional extent of the domain (LENLAT), the zonal "//&
                 "extent (LENLON), and the number of grid points in each "//&
                 "direction are _not_ independent. In MOM the meridional "//&
                 "extent is determined to fit the zonal extent and the "//&
                 "number of grid points, while grid is perfectly isotropic.", &
                 default=.false.)
  call get_param(param_file, mdl, "EQUATOR_REFERENCE", gp%equator_reference, &
                 "If true, the grid is defined to have the equator at the "//&
                 "nearest q- or h- grid point to (-LOWLAT*NJGLOBAL/LENLAT).", &
                 default=.true.)
  call get_param(param_file, mdl, "LAT_ENHANCE_FACTOR", gp%Lat_enhance_factor, &
                 "The amount by which the meridional resolution is "//&
                 "enhanced within LAT_EQ_ENHANCE of the equator.", &
                 units="nondim", default=1.0)
  call get_param(param_file, mdl, "LAT_EQ_ENHANCE", gp%Lat_eq_enhance, &
                 "The latitude range to the north and south of the equator "//&
                 "over which the resolution is enhanced.", units="degrees", &
                 default=0.0)
 
  !   With an isotropic grid, the north-south extent of the domain,
  ! the east-west extent, and the number of grid points in each
  ! direction are _not_ independent.  Here the north-south extent
  ! will be determined to fit the east-west extent and the number of
  ! grid points.  The grid is perfectly isotropic.
  if (gp%equator_reference) then
    ! With the following expression, the equator will always be placed
    ! on either h or q points, in a position consistent with the ratio
    ! GP%south_lat to GP%len_lat.
    jref =  (g%jsg-1) + 0.5*floor(gp%njglobal*((-1.0*gp%south_lat*2.0)/gp%len_lat)+0.5)
    fnref = int_dj_dy(0.0, gp)
  else
    ! The following line sets the reference latitude GP%south_lat at j=js-1 (or -2?)
    jref = (g%jsg-1)
    fnref = int_dj_dy((gp%south_lat*pi/180.0), gp)
  endif
 
  ! These calculations no longer depend on the the order in which they
  ! are performed because they all use the same (poor) starting guess and
  ! iterate to convergence.
  ! Note that the dynamic grid always uses symmetric memory for the global
  ! arrays G%gridLatB and G%gridLonB.
  do j=g%jsg-1,g%jeg
    jd = fnref + (j - jref)
    y_q = find_root(int_dj_dy, dy_dj, gp, jd, 0.0, -1.0*pi_2, pi_2, itt2)
    g%gridLatB(j) = y_q*180.0/pi
    ! if (is_root_pe()) &
    !   write(*, '("J, y_q = ",I4,ES14.4," itts = ",I4)')  j, y_q, itt2
  enddo
  do j=g%jsg,g%jeg
    jd = fnref + (j - jref) - 0.5
    y_h = find_root(int_dj_dy, dy_dj, gp, jd, 0.0, -1.0*pi_2, pi_2, itt1)
    g%gridLatT(j) = y_h*180.0/pi
    ! if (is_root_pe()) &
    !   write(*, '("j, y_h = ",I4,ES14.4," itts = ",I4)')  j, y_h, itt1
  enddo
  do j=jsdb+j_off,jedb+j_off
    jd = fnref + (j - jref)
    y_q = find_root(int_dj_dy, dy_dj, gp, jd, 0.0, -1.0*pi_2, pi_2, itt2)
    do i=isdb,iedb ; yq(i,j-j_off) = y_q ; enddo
    do i=isd,ied ; yv(i,j-j_off) = y_q ; enddo
  enddo
  do j=jsd+j_off,jed+j_off
    jd = fnref + (j - jref) - 0.5
    y_h = find_root(int_dj_dy, dy_dj, gp, jd, 0.0, -1.0*pi_2, pi_2, itt1)
    if ((j >= jsd+j_off) .and. (j <= jed+j_off)) then
      do i=isd,ied ; yh(i,j-j_off) = y_h ; enddo
      do i=isdb,iedb ; yu(i,j-j_off) = y_h ; enddo
    endif
  enddo
 
  ! Determine the longitudes of the various points.
 
  ! These two lines place the western edge of the domain at GP%west_lon.
  iref = (g%isg-1) + gp%niglobal
  fnref = int_di_dx(((gp%west_lon+gp%len_lon)*pi/180.0), gp)
 
  ! These calculations no longer depend on the the order in which they
  ! are performed because they all use the same (poor) starting guess and
  ! iterate to convergence.
  do i=g%isg-1,g%ieg
    id = fnref + (i - iref)
    x_q = find_root(int_di_dx, dx_di, gp, id, 0.0, -4.0*pi, 4.0*pi, itt2)
    g%gridLonB(i) = x_q*180.0/pi
  enddo
  do i=g%isg,g%ieg
    id = fnref + (i - iref) - 0.5
    x_h = find_root(int_di_dx, dx_di, gp, id, 0.0, -4.0*pi, 4.0*pi, itt1)
    g%gridLonT(i) = x_h*180.0/pi
  enddo
  do i=isdb+i_off,iedb+i_off
    id = fnref + (i - iref)
    x_q = find_root(int_di_dx, dx_di, gp, id, 0.0, -4.0*pi, 4.0*pi, itt2)
    do j=jsdb,jedb ; xq(i-i_off,j) = x_q ; enddo
    do j=jsd,jed ; xu(i-i_off,j) = x_q ; enddo
  enddo
  do i=isd+i_off,ied+i_off
    id = fnref + (i - iref) - 0.5
    x_h = find_root(int_di_dx, dx_di, gp, id, 0.0, -4.0*pi, 4.0*pi, itt1)
    do j=jsd,jed ; xh(i-i_off,j) = x_h ; enddo
    do j=jsdb,jedb ; xv(i-i_off,j) = x_h ; enddo
  enddo
 
  do j=jsdb,jedb ; do i=isdb,iedb
    g%geoLonBu(i,j) = xq(i,j)*180.0/pi
    g%geoLatBu(i,j) = yq(i,j)*180.0/pi
    g%dxBu(i,j) = m_to_l*ds_di(xq(i,j), yq(i,j), gp)
    g%dyBu(i,j) = m_to_l*ds_dj(xq(i,j), yq(i,j), gp)
 
    g%areaBu(i,j) = g%dxBu(i,j) * g%dyBu(i,j)
    g%IareaBu(i,j) = 1.0 / (g%areaBu(i,j))
  enddo ; enddo
 
  do j=jsd,jed ; do i=isd,ied
    g%geoLonT(i,j) = xh(i,j)*180.0/pi
    g%geoLatT(i,j) = yh(i,j)*180.0/pi
    g%dxT(i,j) = m_to_l*ds_di(xh(i,j), yh(i,j), gp)
    g%dyT(i,j) = m_to_l*ds_dj(xh(i,j), yh(i,j), gp)
 
    g%areaT(i,j) = g%dxT(i,j)*g%dyT(i,j)
    g%IareaT(i,j) = 1.0 / (g%areaT(i,j))
  enddo ; enddo
 
  do j=jsd,jed ; do i=isdb,iedb
    g%geoLonCu(i,j) = xu(i,j)*180.0/pi
    g%geoLatCu(i,j) = yu(i,j)*180.0/pi
    g%dxCu(i,j) = m_to_l*ds_di(xu(i,j), yu(i,j), gp)
    g%dyCu(i,j) = m_to_l*ds_dj(xu(i,j), yu(i,j), gp)
  enddo ; enddo
 
  do j=jsdb,jedb ; do i=isd,ied
    g%geoLonCv(i,j) = xv(i,j)*180.0/pi
    g%geoLatCv(i,j) = yv(i,j)*180.0/pi
    g%dxCv(i,j) = m_to_l*ds_di(xv(i,j), yv(i,j), gp)
    g%dyCv(i,j) = m_to_l*ds_dj(xv(i,j), yv(i,j), gp)
  enddo ; enddo
 
  if (.not.simple_area) then
    do j=jsdb+1,jed ; do i=isdb+1,ied
      g%areaT(i,j) = m_to_l**2*gp%Rad_Earth**2 * &
          (dl(xq(i-1,j-1),xq(i-1,j),yq(i-1,j-1),yq(i-1,j)) + &
          (dl(xq(i-1,j),xq(i,j),yq(i-1,j),yq(i,j)) +          &
          (dl(xq(i,j),xq(i,j-1),yq(i,j),yq(i,j-1)) +          &
           dl(xq(i,j-1),xq(i-1,j-1),yq(i,j-1),yq(i-1,j-1)))))
    enddo ;enddo
    if ((isdb == isd) .or. (jsdb == jsq)) then
      ! Fill in row and column 1 to calculate the area in the southernmost
      ! and westernmost land cells when we are not using symmetric memory.
      ! The pass_var call updates these values if they are not land cells.
      g%areaT(isd+1,jsd) = g%areaT(isd+1,jsd+1)
      do j=jsd,jed ; g%areaT(isd,j) = g%areaT(isd+1,j) ; enddo
      do i=isd,ied ; g%areaT(i,jsd) = g%areaT(i,jsd+1) ; enddo
      ! Now replace the data in the halos, if value values exist.
      call pass_var(g%areaT,g%Domain)
    endif
    do j=jsd,jed ; do i=isd,ied
      g%IareaT(i,j) = 1.0 / (g%areaT(i,j))
    enddo ; enddo
  endif
 
  call calltree_leave("set_grid_metrics_mercator()")
end subroutine set_grid_metrics_mercator
 
 
!> This function returns the grid spacing in the logical x direction.
function ds_di(x, y, GP)
  real, intent(in) :: x  !< The longitude in question
  real, intent(in) :: y  !< The latitude in question
  type(gps), intent(in) :: gp  !< A structure of grid parameters
  real :: ds_di
  ! Local variables
 
  ds_di = gp%Rad_Earth * cos(y) * dx_di(x,gp)
  ! In general, this might be...
  ! ds_di = GP%Rad_Earth * sqrt( cos(y)*cos(y) * dx_di(x,y,GP)*dx_di(x,y,GP) + &
  !                           dy_di(x,y,GP)*dy_di(x,y,GP))
end function ds_di
 
!> This function returns the grid spacing in the logical y direction.
function ds_dj(x, y, GP)
  real, intent(in) :: x  !< The longitude in question
  real, intent(in) :: y  !< The latitude in question
  type(gps), intent(in) :: gp  !< A structure of grid parameters
  ! Local variables
  real :: ds_dj
 
  ds_dj = gp%Rad_Earth * dy_dj(y,gp)
  ! In general, this might be...
  ! ds_dj = GP%Rad_Earth * sqrt( cos(y)*cos(y) * dx_dj(x,y,GP)*dx_dj(x,y,GP) + &
  !                           dy_dj(x,y,GP)*dy_dj(x,y,GP))
end function ds_dj
 
!> This function returns the contribution from the line integral along one of the four sides of a
!! cell face to the area of a cell, assuming that the sides follow a linear path in latitude and
!! longitude (i.e., on a Mercator grid).
function  dl(x1, x2, y1, y2)
  real, intent(in) :: x1 !< Segment starting longitude, in degrees E.
  real, intent(in) :: x2 !< Segment ending longitude, in degrees E.
  real, intent(in) :: y1 !< Segment ending latitude, in degrees N.
  real, intent(in) :: y2 !< Segment ending latitude, in degrees N.
  ! Local variables
  real :: dl
  real :: r, dy
 
  dy = y2 - y1
 
  if (abs(dy) > 2.5e-8) then
    r = ((1.0 - cos(dy))*cos(y1) + sin(dy)*sin(y1)) / dy
  else
    r = (0.5*dy*cos(y1) + sin(y1))
  endif
  dl = r * (x2 - x1)
 
end function  dl
 
!> This subroutine finds and returns the value of y at which the monotonically increasing
!! function fn takes the value fnval, also returning in ittmax the number of iterations of
!! Newton's method that were used to polish the root.
function find_root( fn, dy_df, GP, fnval, y1, ymin, ymax, ittmax)
  real :: find_root !< The value of y where fn(y) = fnval that will be returned
  real,      external    :: fn    !< The external function whose root is being sought
  real,      external    :: dy_df !< The inverse of the derivative of that function
  type(gps), intent(in)  :: gp  !< A structure of grid parameters
  real,      intent(in)  :: fnval !< The value of fn being sought
  real,      intent(in)  :: y1    !< A first guess for y
  real,      intent(in)  :: ymin  !< The minimum permitted value of y
  real,      intent(in)  :: ymax  !< The maximum permitted value of y
  integer,   intent(out) :: ittmax !< The number of iterations used to polish the root
  ! Local variables
  real :: y, y_next
  real :: ybot, ytop, fnbot, fntop
  integer :: itt
  character(len=256) :: warnmesg
 
  real :: dy_dfn, dy, fny
 
!  Bracket the root.  Do not use the bounding values because the value at the
! function at the bounds could be infinite, as is the case for the Mercator
! grid recursion relation. (I.e., this is a search on an open interval.)
  ybot = y1
  fnbot = fn(ybot,gp) - fnval ; itt = 0
  do while (fnbot > 0.0)
    if ((ybot - 2.0*dy_df(ybot,gp)) < (0.5*(ybot+ymin))) then
      ! Go twice as far as the secant method would normally go.
      ybot = ybot - 2.0*dy_df(ybot,gp)
    else  ! But stay within the open interval!
      ybot = 0.5*(ybot+ymin) ; itt = itt + 1
    endif
    fnbot = fn(ybot,gp) - fnval
 
    if ((itt > 50) .and. (fnbot > 0.0)) then
      write(warnmesg, '("PE ",I2," unable to find bottom bound for grid function. &
        &x = ",ES10.4,", xmax = ",ES10.4,", fn = ",ES10.4,", dfn_dx = ",ES10.4,&
        &", seeking fn = ",ES10.4," - fn = ",ES10.4,".")') &
          pe_here(),ybot,ymin,fn(ybot,gp),dy_df(ybot,gp),fnval, fnbot
      call mom_error(fatal,warnmesg)
    endif
  enddo
 
  ytop = y1
  fntop = fn(ytop,gp) - fnval ; itt = 0
  do while (fntop < 0.0)
    if ((ytop + 2.0*dy_df(ytop,gp)) < (0.5*(ytop+ymax))) then
      ! Go twice as far as the secant method would normally go.
      ytop = ytop + 2.0*dy_df(ytop,gp)
    else ! But stay within the open interval!
      ytop = 0.5*(ytop+ymax) ; itt = itt + 1
    endif
    fntop = fn(ytop,gp) - fnval
 
    if ((itt > 50) .and. (fntop < 0.0)) then
      write(warnmesg, '("PE ",I2," unable to find top bound for grid function. &
        &x = ",ES10.4,", xmax = ",ES10.4,", fn = ",ES10.4,", dfn_dx = ",ES10.4, &
        &", seeking fn = ",ES10.4," - fn = ",ES10.4,".")') &
          pe_here(),ytop,ymax,fn(ytop,gp),dy_df(ytop,gp),fnval,fntop
      call mom_error(fatal,warnmesg)
    endif
  enddo
 
  ! Find the root using a bracketed variant of Newton's method, starting
  ! with a false-positon method first guess.
  if ((fntop < 0.0) .or. (fnbot > 0.0) .or. (ytop < ybot)) then
    write(warnmesg, '("PE ",I2," find_root failed to bracket function. y = ",&
              &2ES10.4,", fn = ",2ES10.4,".")') pe_here(),ybot,ytop,fnbot,fntop
    call mom_error(fatal, warnmesg)
  endif
 
  if (fntop == 0.0) then ; y = ytop ; fny = fntop
  elseif (fnbot == 0.0) then ; y = ybot ; fny = fnbot
  else
    y = (ybot*fntop - ytop*fnbot) / (fntop - fnbot)
    fny = fn(y,gp) - fnval
    if (fny < 0.0) then ; fnbot = fny ; ybot = y
    else ; fntop = fny ; ytop = y ; endif
  endif
 
  do itt=1,50
    dy_dfn = dy_df(y,gp)
 
    dy = -1.0* fny * dy_dfn
    y_next = y + dy
    if ((y_next >= ytop) .or. (y_next <= ybot)) then
      ! The Newton's method estimate has escaped bracketing, so use the
      ! false-position method instead.  The complicated test is to properly
      ! handle the case where the iteration is down to roundoff level differences.
      y_next = y
      if (abs(fntop - fnbot) > epsilon(y) * (abs(fntop) + abs(fnbot))) &
        y_next = (ybot*fntop - ytop*fnbot) / (fntop - fnbot)
    endif
 
    dy = y_next - y
    if (abs(dy) < (2.0*epsilon(y)*(abs(y) + abs(y_next)) + 1.0e-20)) then
      y = y_next ; exit
    endif
    y = y_next
 
    fny = fn(y,gp) - fnval
    if (fny > 0.0) then ; ytop = y ; fntop = fny
    elseif (fny < 0.0) then ; ybot = y ; fnbot = fny
    else ; exit ; endif
 
  enddo
  if (abs(y) < 1e-12) y = 0.0
 
  ittmax = itt
  find_root = y
end function find_root
 
!> This function calculates and returns the value of dx/di, where x is the
!! longitude in Radians, and i is the integral north-south grid index.
function dx_di(x, GP)
  real, intent(in) :: x  !< The longitude in question
  type(gps), intent(in) :: gp  !< A structure of grid parameters
  real :: dx_di
 
  dx_di = (gp%len_lon * 4.0*atan(1.0)) / (180.0 * gp%niglobal)
 
end function dx_di
 
!> This function calculates and returns the integral of the inverse
!! of dx/di to the point x, in radians.
function int_di_dx(x, GP)
  real, intent(in) :: x  !< The longitude in question
  type(gps), intent(in) :: gp  !< A structure of grid parameters
  real :: int_di_dx
 
  int_di_dx = x * ((180.0 * gp%niglobal) / (gp%len_lon * 4.0*atan(1.0)))
 
end function int_di_dx
 
!> This subroutine calculates and returns the value of dy/dj, where y is the
!! latitude in Radians, and j is the integral north-south grid index.
function dy_dj(y, GP)
  real, intent(in) :: y  !< The latitude in question
  type(gps), intent(in) :: gp  !< A structure of grid parameters
  real :: dy_dj
  ! Local variables
  real :: pi            ! 3.1415926... calculated as 4*atan(1)
  real :: c0            ! The constant that converts the nominal y-spacing in
                        ! gridpoints to the nominal spacing in Radians.
  real :: y_eq_enhance  ! The latitude in radians within which the resolution
                        ! is enhanced.
  pi = 4.0*atan(1.0)
  if (gp%isotropic) then
    c0 = (gp%len_lon * pi) / (180.0 * gp%niglobal)
    y_eq_enhance = pi*abs(gp%lat_eq_enhance)/180.0
    if (abs(y) < y_eq_enhance) then
      dy_dj = c0 * (cos(y) / (1.0 + 0.5*cos(y) * (gp%lat_enhance_factor - 1.0) * &
                         (1.0+cos(pi*y/y_eq_enhance)) ))
    else
      dy_dj = c0 * cos(y)
    endif
  else
    c0 = (gp%len_lat * pi) / (180.0 * gp%njglobal)
    dy_dj = c0
  endif
 
end function dy_dj
 
!> This subroutine calculates and returns the integral of the inverse
!! of dy/dj to the point y, in radians.
function int_dj_dy(y, GP)
  real, intent(in) :: y  !< The latitude in question
  type(gps), intent(in) :: gp  !< A structure of grid parameters
  real :: int_dj_dy
  ! Local variables
  real :: i_c0 = 0.0       !   The inverse of the constant that converts the
                           ! nominal spacing in gridpoints to the nominal
                           ! spacing in Radians.
  real :: pi               ! 3.1415926... calculated as 4*atan(1)
  real :: y_eq_enhance     ! The latitude in radians from
                           ! from the equator within which the
                           ! meridional grid spacing is enhanced by
                           ! a factor of GP%lat_enhance_factor.
  real :: r
 
  pi = 4.0*atan(1.0)
  if (gp%isotropic) then
    i_c0 = (180.0 * gp%niglobal) / (gp%len_lon * pi)
    y_eq_enhance = pi*abs(gp%lat_eq_enhance)/180.0
 
    if (y >= 0.0) then
      r = i_c0 * log((1.0 + sin(y))/cos(y))
    else
      r = -1.0 * i_c0 * log((1.0 - sin(y))/cos(y))
    endif
 
    if (y >= y_eq_enhance) then
      r = r + i_c0*0.5*(gp%lat_enhance_factor - 1.0)*y_eq_enhance
    elseif (y <= -y_eq_enhance) then
      r = r - i_c0*0.5*(gp%lat_enhance_factor - 1.0)*y_eq_enhance
    else
      r = r + i_c0*0.5*(gp%lat_enhance_factor - 1.0) * &
              (y + (y_eq_enhance/pi)*sin(pi*y/y_eq_enhance))
    endif
  else
    i_c0 = (180.0 * gp%njglobal) / (gp%len_lat * pi)
    r = i_c0 * y
  endif
 
  int_dj_dy = r
end function int_dj_dy
 
!> Extrapolates missing metric data into all the halo regions.
subroutine extrapolate_metric(var, jh, missing)
  real, dimension(:,:), intent(inout) :: var     !< The array in which to fill in halos
  integer,              intent(in)    :: jh      !< The size of the halos to be filled
  real,       optional, intent(in)    :: missing !< The missing data fill value, 0 by default.
  ! Local variables
  real :: badval
  integer :: i,j
 
  badval = 0.0 ; if (present(missing)) badval = missing
 
  ! Fill in southern halo by extrapolating from the computational domain
  do j=lbound(var,2)+jh,lbound(var,2),-1 ; do i=lbound(var,1),ubound(var,1)
    if (var(i,j)==badval) var(i,j) = 2.0*var(i,j+1)-var(i,j+2)
  enddo ; enddo
 
  ! Fill in northern halo by extrapolating from the computational domain
  do j=ubound(var,2)-jh,ubound(var,2) ; do i=lbound(var,1),ubound(var,1)
    if (var(i,j)==badval) var(i,j) = 2.0*var(i,j-1)-var(i,j-2)
  enddo ; enddo
 
  ! Fill in western halo by extrapolating from the computational domain
  do j=lbound(var,2),ubound(var,2) ; do i=lbound(var,1)+jh,lbound(var,1),-1
    if (var(i,j)==badval) var(i,j) = 2.0*var(i+1,j)-var(i+2,j)
  enddo ; enddo
 
  ! Fill in eastern halo by extrapolating from the computational domain
  do j=lbound(var,2),ubound(var,2) ; do i=ubound(var,1)-jh,ubound(var,1)
    if (var(i,j)==badval) var(i,j) = 2.0*var(i-1,j)-var(i-2,j)
  enddo ; enddo
 
end subroutine extrapolate_metric
 
!> This function implements Adcroft's rule for reciprocals, namely that
!!   Adcroft_Inv(x) = 1/x for |x|>0 or 0 for x=0.
function adcroft_reciprocal(val) result(I_val)
  real, intent(in) :: val  !< The value being inverted.
  real :: i_val            !< The Adcroft reciprocal of val.
 
  i_val = 0.0
  if (val /= 0.0) i_val = 1.0/val
end function adcroft_reciprocal
 
!> Initializes the grid masks and any metrics that come with masks already applied.
!!
!!    Initialize_masks sets mask2dT, mask2dCu, mask2dCv, and mask2dBu to mask out
!! flow over any points which are shallower than Dmin and permit an
!! appropriate treatment of the boundary conditions.  mask2dCu and mask2dCv
!! are 0.0 at any points adjacent to a land point.  mask2dBu is 0.0 at
!! any land or boundary point.  For points in the interior, mask2dCu,
!! mask2dCv, and mask2dBu are all 1.0.
subroutine initialize_masks(G, PF, US)
  type(dyn_horgrid_type),          intent(inout) :: g  !< The dynamic horizontal grid type
  type(param_file_type),           intent(in)    :: pf !< Parameter file structure
  type(unit_scale_type), optional, intent(in)    :: us !< A dimensional unit scaling type
  ! Local variables
  real :: m_to_z_scale ! A unit conversion factor from m to Z.
  real :: m_to_l  ! A unit conversion factor [L m-1 ~> nondim]
  real :: dmin       ! The depth for masking in the same units as G%bathyT [Z ~> m].
  real :: min_depth  ! The minimum ocean depth in the same units as G%bathyT [Z ~> m].
  real :: mask_depth ! The depth shallower than which to mask a point as land [Z ~> m].
  character(len=40)  :: mdl = "MOM_grid_init initialize_masks"
  integer :: i, j
 
  call calltree_enter("initialize_masks(), MOM_grid_initialize.F90")
  m_to_z_scale = 1.0 ; if (present(us)) m_to_z_scale = us%m_to_Z
  m_to_l = 1.0 ; if (present(us)) m_to_l = us%m_to_L
 
  call get_param(pf, 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_scale)
  call get_param(pf, mdl, "MASKING_DEPTH", mask_depth, &
                 "The depth below which to mask points as land points, for which all "//&
                 "fluxes are zeroed out. MASKING_DEPTH is ignored if negative.", &
                 units="m", default=-9999.0, scale=m_to_z_scale)
 
  dmin = min_depth
  if (mask_depth>=0.) dmin = mask_depth
 
  g%mask2dCu(:,:) = 0.0 ; g%mask2dCv(:,:) = 0.0 ; g%mask2dBu(:,:) = 0.0
 
  ! Construct the h-point or T-point mask
  do j=g%jsd,g%jed ; do i=g%isd,g%ied
    if (g%bathyT(i,j) <= dmin) then
      g%mask2dT(i,j) = 0.0
    else
      g%mask2dT(i,j) = 1.0
    endif
  enddo ; enddo
 
  do j=g%jsd,g%jed ; do i=g%isd,g%ied-1
    if ((g%bathyT(i,j) <= dmin) .or. (g%bathyT(i+1,j) <= dmin)) then
      g%mask2dCu(i,j) = 0.0
    else
      g%mask2dCu(i,j) = 1.0
    endif
  enddo ; enddo
 
  do j=g%jsd,g%jed-1 ; do i=g%isd,g%ied
    if ((g%bathyT(i,j) <= dmin) .or. (g%bathyT(i,j+1) <= dmin)) then
      g%mask2dCv(i,j) = 0.0
    else
      g%mask2dCv(i,j) = 1.0
    endif
  enddo ; enddo
 
  do j=g%jsd,g%jed-1 ; do i=g%isd,g%ied-1
    if ((g%bathyT(i+1,j) <= dmin) .or. (g%bathyT(i+1,j+1) <= dmin) .or. &
        (g%bathyT(i,j) <= dmin) .or. (g%bathyT(i,j+1) <= dmin)) then
      g%mask2dBu(i,j) = 0.0
    else
      g%mask2dBu(i,j) = 1.0
    endif
  enddo ; enddo
 
  call pass_var(g%mask2dBu, g%Domain, position=corner)
  call pass_vector(g%mask2dCu, g%mask2dCv, g%Domain, to_all+scalar_pair, cgrid_ne)
 
  do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB
    g%dy_Cu(i,j) = g%mask2dCu(i,j) * g%dyCu(i,j)
    g%areaCu(i,j) = g%dxCu(i,j) * g%dy_Cu(i,j)
    g%IareaCu(i,j) = g%mask2dCu(i,j) * adcroft_reciprocal(g%areaCu(i,j))
  enddo ; enddo
 
  do j=g%JsdB,g%JedB ; do i=g%isd,g%ied
    g%dx_Cv(i,j) = g%mask2dCv(i,j) * g%dxCv(i,j)
    g%areaCv(i,j) = g%dyCv(i,j) * g%dx_Cv(i,j)
    g%IareaCv(i,j) = g%mask2dCv(i,j) * adcroft_reciprocal(g%areaCv(i,j))
  enddo ; enddo
 
  call calltree_leave("initialize_masks()")
end subroutine initialize_masks
 
!> \namespace mom_grid_initialize
!!
!!  The metric terms have the form Dzp, IDzp, or DXDYp, where z can
!!  be X or Y, and p can be q, u, v, or h.  z describes the direction
!!  of the metric, while p describes the location.  IDzp is the
!!  inverse of Dzp, while DXDYp is the product of DXp and DYp except
!!  that areaT is calculated analytically from the latitudes and
!!  longitudes of the surrounding q points.
!!
!!    On a sphere, a variety of grids can be implemented by defining
!!  analytic expressions for dx_di, dy_dj (where x and y are latitude
!!  and longitude, and i and j are grid indices) and the expressions
!!  for the integrals of their inverses in the four subroutines
!!  dy_dj, Int_dj_dy, dx_di, and Int_di_dx.
!!
!!    initialize_masks sets up land masks based on the depth field.
!!  The one argument is the minimum ocean depth.  Depths that are
!!  less than this are interpreted as land points.
 
end module mom_grid_initialize