MOM_grid_initialize.F90

Go to the documentation of this file.

module mom_grid_initialize
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg, November 1998 - June 2002                      *
!*                                                                     *
!*    This program contains 2 externally callable subroutines.         *
!*  set_grid_metrics calculates the various metric terms that are used *
!*  by MOM.  This routine is intended to be modified by the user to    *
!*  enable the use of any general orthogonal grid.  initialize_masks   *
!*  initializes the land masks; it is in this file because it a key    *
!*  part of the physical grid description.                             *
!*                                                                     *
!*    This subroutine is also used by MOM-related preprocessing and    *
!*  postprocessing codes.                                              *
!*                                                                     *
!*    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.                     *
!*                                                                     *
!*    Macros written all in capital letters are from MOM_memory.h.     *
!*                                                                     *
!*     A small fragment of the C-grid is shown below:                  *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q, dxBu, IdxBu, dyBu, IdyBu, etc.        *
!*    j+1  > o > o >   At ^:  v, dxCv, IdxCv, dyCv, IdyCv, etc.        *
!*    j    x ^ x ^ x   At >:  u, dxCu, IdxCu, dyCu, IdyCu, etc.        *
!*    j    > o > o >   At o:  h, dxT, IdxT, dyT, IdyT, areaT, etc.     *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

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 : read_data, slasher, file_exists
use mom_io, only : corner, north_face, east_face

use mpp_domains_mod, only : mpp_get_domain_extents, mpp_deallocate_domain

implicit none ; private

public set_grid_metrics, initialize_masks, adcroft_reciprocal

type, public :: gps ; private
  real :: len_lon
  real :: len_lat
  real :: west_lon
  real :: south_lat
  real :: rad_earth
  real :: lat_enhance_factor
  real :: lat_eq_enhance
  logical :: isotropic
  logical :: equator_reference
  integer :: niglobal, njglobal         ! Duplicates of niglobal and njglobal from MOM_dom
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)
  type(dyn_horgrid_type), intent(inout) :: G          !< The dynamic horizontal grid type
  type(param_file_type), intent(in)    :: param_file  !< Parameter file structure
! Arguments:
!  (inout)   G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.

!    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.
! 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 \n"//&
                 "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.)

  ! These are defaults that may be changed in the next select block.
  g%x_axis_units = "degrees_E" ; g%y_axis_units = "degrees_N"
  select case (trim(config))
    case ("mosaic");    call set_grid_metrics_from_mosaic(g, param_file)
    case ("cartesian"); call set_grid_metrics_cartesian(g, param_file)
    case ("spherical"); call set_grid_metrics_spherical(g, param_file)
    case ("mercator");  call set_grid_metrics_mercator(g, param_file)
    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)
  call calltree_leave("set_derived_metrics()")

  if (debug) call grid_metrics_chksum('MOM_grid_init/set_grid_metrics',g)

  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)
  character(len=*),      intent(in) :: parent  !< A string identifying the caller
  type(dyn_horgrid_type), intent(in) :: G      !< The dynamic horizontal grid type

  integer :: halo

  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)

  call uvchksum(trim(parent)//': dxC[uv]', g%dxCu, g%dyCv, g%HI, haloshift=halo)

  call uvchksum(trim(parent)//': dxC[uv]', &
                g%dyCu, g%dxCv, g%HI, haloshift=halo)

  call bchksum_pair(trim(parent)//': dxB[uv]', &
                    g%dxBu, g%dyBu, g%HI, haloshift=halo)

  call hchksum_pair(trim(parent)//': Id[xy]T', &
                    g%IdxT, g%IdyT, g%HI, haloshift=halo)

  call uvchksum(trim(parent)//': Id[xy]C[uv]', &
                g%IdxCu, g%IdyCv, g%HI, haloshift=halo)

  call uvchksum(trim(parent)//': Id[xy]C[uv]', &
                g%IdyCu, g%IdxCv, g%HI, haloshift=halo)

  call bchksum_pair(trim(parent)//': Id[xy]B[uv]', &
                    g%IdxBu, g%IdyBu, g%HI, haloshift=halo)

  call hchksum(g%areaT, trim(parent)//': areaT',g%HI, haloshift=halo)
  call bchksum(g%areaBu, trim(parent)//': areaBu',g%HI, haloshift=halo)

  call hchksum(g%IareaT, trim(parent)//': IareaT',g%HI, haloshift=halo)
  call bchksum(g%IareaBu, trim(parent)//': IareaBu',g%HI, haloshift=halo)

  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

! ------------------------------------------------------------------------------

!>  set_grid_metrics_from_mosaic sets the grid metrics from a mosaic file.
subroutine set_grid_metrics_from_mosaic(G, param_file)
  type(dyn_horgrid_type), intent(inout) :: G           !< The dynamic horizontal grid type
  type(param_file_type), intent(in)     :: param_file  !< Parameter file structure
!   This subroutine sets the grid metrics from a mosaic file.
!
! Arguments:
!  (inout)   G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.

  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
  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
  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")

  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, "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%use_io_layout = g%domain%use_io_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

  if (sgdom%use_io_layout) &
    call mom_define_io_domain(sgdom%mpp_domain, sgdom%io_layout)
  deallocate(exni)
  deallocate(exnj)

! Read X from the supergrid
  tmpz(:,:) = 999.
  call read_data(filename, 'x', tmpz, domain=sgdom%mpp_domain, 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%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 read_data(filename, 'y', tmpz, domain=sgdom%mpp_domain, 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 read_data(filename,'dx',tmpv,domain=sgdom%mpp_domain,position=north_face)
  call read_data(filename,'dy',tmpu,domain=sgdom%mpp_domain,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 read_data(filename, 'area', tmpt, domain=sgdom%mpp_domain)
  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) = dxt(i,j) ; g%dyT(i,j) = dyt(i,j) ; g%areaT(i,j) = areat(i,j)
  enddo ; enddo
  do i=g%IsdB,g%IedB ; do j=g%jsd,g%jed
    g%dxCu(i,j) = dxcu(i,j) ; g%dyCu(i,j) = dycu(i,j)
  enddo ; enddo
  do i=g%isd,g%ied ; do j=g%JsdB,g%JedB
    g%dxCv(i,j) = dxcv(i,j) ; g%dyCv(i,j) = dycv(i,j)
  enddo ; enddo
  do i=g%IsdB,g%IedB ; do j=g%JsdB,g%JedB
    g%dxBu(i,j) = dxbu(i,j) ; g%dyBu(i,j) = dybu(i,j) ; g%areaBu(i,j) = 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


! ------------------------------------------------------------------------------

subroutine set_grid_metrics_cartesian(G, param_file)
  type(dyn_horgrid_type), intent(inout) :: G           !< The dynamic horizontal grid type
  type(param_file_type), intent(in)     :: param_file  !< Parameter file structure

! Arguments:
!  (inout)   G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.

!    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.
  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 in m.
  real :: I_dx, I_dy                   ! Inverse grid spacings in m.
  real :: PI
  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")

  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 \n"//&
                 "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 \n"//&
                 "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))
  else if (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) = dx_everywhere ; g%IdxBu(i,j) = i_dx
    g%dyBu(i,j) = dy_everywhere ; g%IdyBu(i,j) = i_dy
    g%areaBu(i,j) = dx_everywhere * dy_everywhere ; g%IareaBu(i,j) = 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) = dx_everywhere ; g%IdxT(i,j) = i_dx
    g%dyT(i,j) = dy_everywhere ; g%IdyT(i,j) = i_dy
    g%areaT(i,j) = dx_everywhere * dy_everywhere ; g%IareaT(i,j) = 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) = dx_everywhere ; g%IdxCu(i,j) = i_dx
    g%dyCu(i,j) = dy_everywhere ; g%IdyCu(i,j) = 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) = dx_everywhere ; g%IdxCv(i,j) = i_dx
    g%dyCv(i,j) = dy_everywhere ; g%IdyCv(i,j) = i_dy
  enddo ; enddo

  call calltree_leave("set_grid_metrics_cartesian()")
end subroutine set_grid_metrics_cartesian

! ------------------------------------------------------------------------------

subroutine set_grid_metrics_spherical(G, param_file)
  type(dyn_horgrid_type), intent(inout) :: G           !< The dynamic horizontal grid type
  type(param_file_type), intent(in)     :: param_file  !< Parameter file structure

! Arguments:
!  (inout)   G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.

!    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.
  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
  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")

!    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) = g%Rad_Earth * cos( g%geoLatBu(i,j)*pi_180 ) * dl_di
!   G%dxBu(I,J) = G%Rad_Earth * dLon*PI_180 * COS( G%geoLatBu(I,J)*PI_180 )
    g%dyBu(i,j) = 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) = g%Rad_Earth * cos( g%geoLatCv(i,j)*pi_180 ) * dl_di
!   G%dxCv(i,J) = G%Rad_Earth * (dLon*PI_180) * COS( G%geoLatCv(i,J)*PI_180 )
    g%dyCv(i,j) = 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) = g%Rad_Earth * cos( g%geoLatCu(i,j)*pi_180 ) * dl_di
!   G%dxCu(I,j) = G%Rad_Earth * dLon*PI_180 * COS( latitude )
    g%dyCu(i,j) = 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) = 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) = 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) = 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

! ------------------------------------------------------------------------------

subroutine set_grid_metrics_mercator(G, param_file)
  type(dyn_horgrid_type), intent(inout) :: G           !< The dynamic horizontal grid type
  type(param_file_type), intent(in)     :: param_file  !< Parameter file structure

! Arguments:
!  (inout)   G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.

!    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.
  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


!   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.

!  The metric terms within the computational domain are set here.
  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.
  integer :: itt1, itt2
  logical :: debug = .false., simple_area = .true.
  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, IsdB, IedB, JsdB, JedB

  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")

!    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 \n"//&
                 "a Mercator grid) is used. With an isotropic grid, the \n"//&
                 "meridional extent of the domain (LENLAT), the zonal \n"//&
                 "extent (LENLON), and the number of grid points in each \n"//&
                 "direction are _not_ independent. In MOM the meridional \n"//&
                 "extent is determined to fit the zonal extent and the \n"//&
                 "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 \n"//&
                 "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 \n"//&
                 "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 \n"//&
                 "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) = ds_di(xq(i,j), yq(i,j), gp)
    g%dyBu(i,j) = 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) = ds_di(xh(i,j), yh(i,j), gp)
    g%dyT(i,j) = 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) = ds_di(xu(i,j), yu(i,j), gp)
    g%dyCu(i,j) = 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) = ds_di(xv(i,j), yv(i,j), gp)
    g%dyCv(i,j) = 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) = 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


function ds_di(x, y, GP)
  real, intent(in) :: x, y
  type(gps), intent(in) :: GP
  real :: ds_di
! This function returns the grid spacing in the logical x direction.
! Arguments: x - The latitude in question.
!  (in)      y - The longitude in question.
  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

function ds_dj(x, y, GP)
  real, intent(in) :: x, y
  type(gps), intent(in) :: GP
  real :: ds_dj
! This function returns the grid spacing in the logical y direction.
! Arguments: x - The latitude in question.
!  (in)      y - The longitude in question.
  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


function  dl(x1, x2, y1, y2)
  real, intent(in) :: x1, x2, y1, y2
  real :: dL
!  This subroutine calculates 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 long-
! itude (i.e., on a Mercator grid).
! Argumnts: x1 - Segment starting longitude.
!  (in)     x2 - Segment ending longitude.
!  (in)     y1 - Segment ending latitude.
!  (in)     y2 - Segment ending latitude.
  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

function find_root( fn, dy_df, GP, fnval, y1, ymin, ymax, ittmax)
  real :: find_root
  real, external :: fn, dy_df
  type(gps), intent(in) :: GP
  real, intent(in) :: fnval, y1, ymin, ymax
  integer, intent(out) :: ittmax
  real :: y, y_next
! 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.
  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

function dx_di(x, GP)
  real, intent(in) :: x
  type(gps), intent(in) :: GP
  real :: dx_di
! This subroutine calculates and returns the value of dx/di, where
! x is the longitude in Radians, and i is the integral north-south
! grid index.

  dx_di = (gp%len_lon * 4.0*atan(1.0)) / (180.0 * gp%niglobal)

end function dx_di

function int_di_dx(x, GP)
  real, intent(in) :: x
  type(gps), intent(in) :: GP
  real :: Int_di_dx
! This subroutine calculates and returns the integral of the inverse
! of dx/di to the point x, in radians.

  int_di_dx = x * ((180.0 * gp%niglobal) / (gp%len_lon * 4.0*atan(1.0)))

end function int_di_dx

function dy_dj(y, GP)
  real, intent(in) :: y
  type(gps), intent(in) :: GP
  real :: dy_dj
! 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.
  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

function int_dj_dy(y, GP)
  real, intent(in) :: y
  type(gps), intent(in) :: GP
  real :: Int_dj_dy
! This subroutine calculates and returns the integral of the inverse
! of dy/dj to the point y, in radians.
  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
    else if (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

! ------------------------------------------------------------------------------

! ------------------------------------------------------------------------------

!> extrapolate_metric 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.
  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

!> initialize_masks initializes the grid masks and any metrics that come
!!    with masks already applied.
subroutine initialize_masks(G, PF)
  type(dyn_horgrid_type), intent(inout) :: G   !< The dynamic horizontal grid type
  type(param_file_type), intent(in)     :: PF  !< Parameter file structure

!    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.

  real :: Dmin, min_depth, mask_depth
  character(len=40)  :: mdl = "MOM_grid_init initialize_masks"
  integer :: i, j

  call calltree_enter("initialize_masks(), MOM_grid_initialize.F90")
  call get_param(pf, mdl, "MINIMUM_DEPTH", min_depth, &
                 "If MASKING_DEPTH is unspecified, then anything shallower than\n"//&
                 "MINIMUM_DEPTH is assumed to be land and all fluxes are masked out.\n"//&
                 "If MASKING_DEPTH is specified, then all depths shallower than\n"//&
                 "MINIMUM_DEPTH but deeper than MASKING_DEPTH are rounded to MINIMUM_DEPTH.", &
                 units="m", default=0.0)
  call get_param(pf, mdl, "MASKING_DEPTH", mask_depth, &
                 "The depth below which to mask points as land points, for which all\n"//&
                 "fluxes are zeroed out. MASKING_DEPTH is ignored if negative.", &
                 units="m", default=-9999.0)

  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%dy_Cu_obc(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%dx_Cv_obc(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

end module mom_grid_initialize