!--- cpp defs common to all files
#include "cpl_cppdefs"

module ccc_esmf
  !***************************************************************
  !--- CCCma specific routines that use ESMF procedures
  !---
  !--- Larry Solheim  ...Sep,2015
  !***************************************************************

  use ESMF

  use com_cpl
  use grid_comm
  use strings, only: lowerc, upperc, split_on_delim
  use ccc_netcdf
  use ccc_io
  use store

  implicit none

  public

  !--- Generic interface for the various flavors of regrid_Field
  interface regrid_Field
    module procedure regrid_FieldById
    module procedure regrid_FieldByName
  end interface

  !--- Generic interface for a routine that will extract grid coordinate
  !--- data from a grid that is defined in grid_list(:)
  interface get_grid_coord
    module procedure get_grid_coord_2di4
    module procedure get_grid_coord_2di8
    module procedure get_grid_coord_2dr4
    module procedure get_grid_coord_2dr8
  end interface
  private :: get_grid_coord_2di4
  private :: get_grid_coord_2di8
  private :: get_grid_coord_2dr4
  private :: get_grid_coord_2dr8

  !--- regrid_field_cv is a wrapper for a conservative fracarea remap
  public :: regrid_field_cv

  !--- Local verbosity flag (a larger value for verbose means more output)
  integer, private :: verbose=1

  !--- Keep track of the number of MPI tasks used by the coupler
  integer(ESMF_KIND_I), save :: cpl_PetCount=1

  !--- Default log files will have names of the form PET0.default_log_file
  character(128) :: default_log_file = " "

  !--- A list of alarms to use with the coupler clock
  !--- Once allocated, this must have the same size as the event array since
  !--- these alarms are to be associated, one to one, with each event
  type(ESMF_Alarm), pointer, save :: alarm(:) => null()

  type cpl_grid_t
    character(128) ::     title = " "
    character(128) :: file_name = " "
    integer        ::        nx = -1
    integer        ::        ny = -1
    integer        ::      size = -1
    integer, dimension(2,2)  :: edges = -1
    integer(ESMF_KIND_I4), pointer, dimension(:,:) :: mask => null()
    type(ESMF_Grid), pointer :: grid => null()
  end type cpl_grid_t

  !--- A list of all grids used by component models
  integer(ESMF_KIND_I), parameter :: max_grid_list = 100
  integer(ESMF_KIND_I), save :: n_grid_list = 0
  type(ESMF_Grid),  pointer, save :: grid_list(:) => null()
  type(cpl_grid_t), pointer, save :: cpl_grid(:) => null()

  !--- The default grid for each component model is assigned an index in this
  !--- list (currently id_atm_grid, id_ocn_grid, id_ice_grid).
  !--- Additional grids may be defined and added to this list as needed.
  integer, save :: id_atm_grid=-1, id_ocn_grid=-1, id_ice_grid=-1, id_canom4_grid=-1

  type regrid_t
    integer ::   regrid_id  = -1
    integer :: src_grid_id  = -1
    integer :: dst_grid_id  = -1
    logical :: use_src_mask = .false.
    logical :: use_dst_mask = .false.
    character(32)  :: regrid_method = " "
    character(32)  :: norm_type = "fracarea"
    character(32)  :: line_type = "spherical"
    character(128) :: src_grid_name = " "
    character(128) :: dst_grid_name = " "
    character(128) :: src_grid_file_name = " "
    character(128) :: dst_grid_file_name = " "
    character(128) :: wgt_grid_file_name = " "
    integer(ESMF_KIND_I), pointer :: src_mask_values(:) => null()
    integer(ESMF_KIND_I), pointer :: dst_mask_values(:) => null()
    real(ESMF_KIND_R8),   pointer :: src_frac(:,:) => null()
    real(ESMF_KIND_R8),   pointer :: dst_frac(:,:) => null()
    real(ESMF_KIND_R8),   pointer :: src_area(:,:) => null()
    real(ESMF_KIND_R8),   pointer :: dst_area(:,:) => null()
    type(ESMF_RouteHandle), pointer :: rh => null()
  end type regrid_t

  integer, save :: n_regrid_list = 0
  integer, parameter :: max_regrid_list = 20
  type(regrid_t), save :: regrid_list(max_regrid_list)
  type(ESMF_RouteHandle), save, target :: regrid_rh_list(max_regrid_list)

  !--- Logical variables passed into ESMF routines need to be of the same kind
  !--- as the default logical variable that was compiled into the ESMF library
  !--- On AIX, with the xlf compiler the default logical kind is the same as
  !--- the default integer kind (ie if ESMF is compiled with -qintsize=4 then
  !--- the default logical variable is also 4 bytes). This may not be true on
  !--- another machine and/or compiler.
  !--- Define ESMF_KIND_L to be the ESMF default logical kind
  integer, private, parameter :: ESMF_KIND_L = ESMF_KIND_I

  !--- Define True and False parameters to use in calls to ESMF routines
  !--- NOTE: These cannot be named esmf_true and esmf_false because these names
  !---       are defined and made public by the ESMF module. Moreover the
  !---       variables ESMF_TRUE and ESMF_FALSE are of type ESMF_Logical which
  !---       is not the same as type logical.
  logical(ESMF_KIND_L), private, parameter :: T_esmf=.true., F_esmf=.false.

  !--- Define some local parameters
  integer(ESMF_KIND_I), private, parameter :: zero=0
  integer(ESMF_KIND_I), private, parameter ::  one=1
  integer(ESMF_KIND_I), private, parameter ::  two=2
  real(ESMF_KIND_R8), private, parameter   :: pi = 3.141592653589793_8
  real(ESMF_KIND_R8), private, parameter   :: rad2deg = 180.0_8/pi

  contains

    !***************************************************************************
    !--- Extract 1D lon coordinate values from the atm grid definition
    !***************************************************************************
    subroutine get_atm_lon_1d(atm_lon, overlap)
      real(kind=8), pointer :: atm_lon(:)
      logical, intent(in), optional :: overlap

      !--- Local
      logical :: add_overlap
      real(kind=8), pointer :: mid_lon(:,:) => null()

      if ( present(overlap) ) then
        add_overlap = overlap
      else
        add_overlap = .false.
      endif

      if (id_atm_grid < 0 ) then
        write(6,*)"get_atm_lon_1d: ATM grid is not defined"
        call xit("get_atm_lon_1d",-1)
      endif

      !--- Retrieve AGCM cell center lon (dim 1) coordinate values
      call get_grid_coord(mid_lon, id_atm_grid, "center", 1)

      if ( associated(atm_lon) ) deallocate(atm_lon)
      if (add_overlap) then
        if ( nlon_a == size(mid_lon,dim=1) ) then
          !--- Add the overlap longitude
          allocate( atm_lon(nlon_a+olap_a) )
          atm_lon(1:nlon_a) = mid_lon(:,1)
          if ( olap_a == 1 ) then
            atm_lon(nlon_a+1) = mid_lon(1,1) + 360.0_8
          else
            atm_lon(nlon_a+1:nlon_a+olap_a) = mid_lon(1:olap_a,1) + 360.0_8
          endif
        else
          write(6,*)"get_atm_lon_1d: Shape mismatch assigning atm_lon with overlap lon"
          write(6,*)"    expected shape = ",nlon_a
          write(6,*)"       found shape = ",size(mid_lon,dim=1)
          call xit("get_atm_lon_1d",-2)
        endif
      else
        allocate( atm_lon(size(mid_lon,dim=1)) )
        atm_lon = mid_lon(:,1)
      endif
      deallocate(mid_lon)

    end subroutine get_atm_lon_1d

    !***************************************************************************
    !--- Extract 1D lat coordinate values from the atm grid definition
    !***************************************************************************
    subroutine get_atm_lat_1d(atm_lat)
      real(kind=8), pointer, intent(inout) :: atm_lat(:)

      !--- Local
      real(kind=8), pointer :: mid_lat(:,:) => null()

      if (id_atm_grid < 0 ) then
        write(6,*)"get_atm_lat: ATM grid is not defined"
        call xit("get_atm_lat_1d",-1)
      endif

      !--- Retrieve AGCM cell center lat (dim 2) coordinate values
      call get_grid_coord(mid_lat, id_atm_grid, "center", 2)

      if ( associated(atm_lat) ) deallocate(atm_lat)
      allocate( atm_lat(size(mid_lat,dim=2)) )
      atm_lat = mid_lat(1,:)
      deallocate(mid_lat)

    end subroutine get_atm_lat_1d

    !***************************************************************************
    !--- Extract the 2D lon coordinate values from the atm grid definition
    !***************************************************************************
    subroutine get_atm_lon(atm_lon, overlap)
      real(kind=8), pointer :: atm_lon(:,:)
      logical, intent(in), optional :: overlap

      !--- Local
      logical :: add_overlap
      real(kind=8), pointer :: mid_lon(:,:) => null()

      if ( present(overlap) ) then
        add_overlap = overlap
      else
        add_overlap = .false.
      endif

      if (id_atm_grid < 0 ) then
        write(6,*)"get_atm_lon: ATM grid is not defined"
        call xit("get_atm_lon",-1)
      endif

      !--- Retrieve AGCM cell center lon (dim 1) coordinate values
      call get_grid_coord(mid_lon, id_atm_grid, "center", 1)

      if ( associated(atm_lon) ) deallocate(atm_lon)
      if (add_overlap) then
        if ( nlon_a == size(mid_lon,dim=1) ) then
          !--- Add the overlap longitude
          allocate( atm_lon(nlon_a+olap_a,size(mid_lon,dim=2)) )
          atm_lon(1:nlon_a,:) = mid_lon
          if ( olap_a == 1 ) then
            atm_lon(nlon_a+1,:) = mid_lon(1,:) + 360.0_8
          else
            atm_lon(nlon_a+1:nlon_a+olap_a,:) = mid_lon(1:olap_a,:) + 360.0_8
          endif
        else
          write(6,*)"get_atm_lon: Shape mismatch assigning atm_lon with overlap lon"
          write(6,*)"    expected shape = ",nlon_a,nlat_a
          write(6,*)"       found shape = ",shape(mid_lon)
          call xit("get_atm_lon",-2)
        endif
      else
        allocate( atm_lon(size(mid_lon,dim=1),size(mid_lon,dim=2)) )
        atm_lon = mid_lon
      endif
      deallocate(mid_lon)

    end subroutine get_atm_lon

    !***************************************************************************
    !--- Extract the 2D lat coordinate values from the atm grid definition
    !***************************************************************************
    subroutine get_atm_lat(atm_lat, overlap)
      real(kind=8), pointer, intent(inout) :: atm_lat(:,:)
      logical, intent(in), optional :: overlap

      !--- Local
      logical :: add_overlap
      real(kind=8), pointer :: mid_lat(:,:) => null()

      if ( present(overlap) ) then
        add_overlap = overlap
      else
        add_overlap = .false.
      endif

      if (id_atm_grid < 0 ) then
        write(6,*)"get_atm_lat: ATM grid is not defined"
        call xit("get_atm_lat",-1)
      endif

      !--- Retrieve AGCM cell center lat (dim 2) coordinate values
      call get_grid_coord(mid_lat, id_atm_grid, "center", 2)

      if ( associated(atm_lat) ) deallocate(atm_lat)
      if (add_overlap) then
        if ( nlon_a == size(mid_lat,dim=1) ) then
          !--- Add the overlap latgitude
          allocate( atm_lat(nlon_a+olap_a,size(mid_lat,dim=2)) )
          atm_lat(1:nlon_a,:) = mid_lat
          atm_lat(nlon_a+1:nlon_a+olap_a,:) = mid_lat(1:olap_a,:)
        else
          write(6,*)"get_atm_lat: Shape mismatch assigning atm_lat with overlap lat"
          write(6,*)"    expected shape = ",nlon_a,nlat_a
          write(6,*)"       found shape = ",shape(mid_lat)
          call xit("get_atm_lat",-2)
        endif
      else
        allocate( atm_lat(size(mid_lat,dim=1),size(mid_lat,dim=2)) )
        atm_lat = mid_lat
      endif
      deallocate(mid_lat)

    end subroutine get_atm_lat

    !***************************************************************************
    !--- Extract the 2D lon coordinate values from the ocn grid definition
    !***************************************************************************
    subroutine get_ocn_lon(ocn_lon, overlap)
      real(kind=8), pointer :: ocn_lon(:,:)
      logical, intent(in), optional :: overlap

      !--- Local
      logical :: add_overlap
      real(kind=8), pointer :: mid_lon(:,:) => null()

      if ( present(overlap) ) then
        add_overlap = overlap
      else
        add_overlap = .false.
      endif

      if (id_ocn_grid < 0 ) then
        write(6,*)"get_ocn_lon: OCN grid is not defined"
        call xit("get_ocn_lon",-1)
      endif

      !--- Retrieve OCN cell center lon (dim 1) coordinate values
      call get_grid_coord(mid_lon, id_ocn_grid, "center", 1)

      if ( associated(ocn_lon) ) deallocate(ocn_lon)
      if (add_overlap) then
        if ( nlon_o == size(mid_lon,dim=1) ) then
          !--- Add the overlap longitude
          allocate( ocn_lon(nlon_o+2,size(mid_lon,dim=2)) )
          !--- First assign the interior points (excluding overlap lons)
          ocn_lon(2:nlon_o+1,:) = mid_lon
          !--- Then assign the extra longitudes
          ocn_lon(1,:) = mid_lon(nlon_o,:)
          ocn_lon(nlon_o+2,:) = mid_lon(1,:)
        else
          write(6,*)"get_ocn_lon: Shape mismatch assigning ocn_lon with overlap lon"
          write(6,*)"    expected shape = ",nlon_o,nlat_o
          write(6,*)"       found shape = ",shape(mid_lon)
          call xit("get_ocn_lon",-2)
        endif
      else
        allocate( ocn_lon(size(mid_lon,dim=1),size(mid_lon,dim=2)) )
        ocn_lon = mid_lon
      endif
      deallocate(mid_lon)

    end subroutine get_ocn_lon

    !***************************************************************************
    !--- Extract the 2D lat coordinate values from the ocn grid definition
    !***************************************************************************
    subroutine get_ocn_lat(ocn_lat, overlap)
      real(kind=8), pointer, intent(inout) :: ocn_lat(:,:)
      logical, intent(in), optional :: overlap

      !--- Local
      logical :: add_overlap
      real(kind=8), pointer :: mid_lat(:,:) => null()

      if ( present(overlap) ) then
        add_overlap = overlap
      else
        add_overlap = .false.
      endif

      if (id_ocn_grid < 0 ) then
        write(6,*)"get_ocn_lat: OCN grid is not defined"
        call xit("get_ocn_lat",-1)
      endif

      !--- Retrieve OCN cell center lat (dim 2) coordinate values
      call get_grid_coord(mid_lat, id_ocn_grid, "center", 2)

      if ( associated(ocn_lat) ) deallocate(ocn_lat)
      if (add_overlap) then
        if ( nlon_o == size(mid_lat,dim=1) ) then
          !--- Add the overlap longitude
          allocate( ocn_lat(nlon_o+2,size(mid_lat,dim=2)) )
          !--- First assign the interior points (excluding overlap lons)
          ocn_lat(2:nlon_o+1,:) = mid_lat
          !--- Then assign the extra longitudes
          ocn_lat(1,:) = mid_lat(nlon_o,:)
          ocn_lat(nlon_o+2,:) = mid_lat(1,:)
        else
          write(6,*)"get_ocn_lat: Shape mismatch assigning ocn_lat with overlap lon"
          write(6,*)"    expected shape = ",nlon_o,nlat_o
          write(6,*)"       found shape = ",shape(mid_lat)
          call xit("get_ocn_lat",-2)
        endif
      else
        allocate( ocn_lat(size(mid_lat,dim=1),size(mid_lat,dim=2)) )
        ocn_lat = mid_lat
      endif
      deallocate(mid_lat)

    end subroutine get_ocn_lat

    !***************************************************************************
    !--- Extract the 2D mask values from the ocn grid definition
    !***************************************************************************
    subroutine get_ocn_mask(ocn_mask, overlap)
      integer, pointer, intent(inout) :: ocn_mask(:,:)
      logical, intent(in), optional :: overlap

      !--- Local
      logical :: add_overlap
      real(kind=8), pointer :: tmp_mask(:,:) => null()

      if ( present(overlap) ) then
        add_overlap = overlap
      else
        add_overlap = .false.
      endif

      if (id_ocn_grid < 0 ) then
        write(6,*)"get_ocn_mask: OCN grid is not defined"
        call xit("get_ocn_mask",-1)
      endif

      !--- Retrieve OCN cell center mask values
      call get_grid_item(tmp_mask, id_ocn_grid, "center", "mask")

      if ( associated(ocn_mask) ) deallocate(ocn_mask)
      if (add_overlap) then
        if ( nlon_o == size(tmp_mask,dim=1) ) then
          !--- Add the overlap longitude
          allocate( ocn_mask(nlon_o+2,size(tmp_mask,dim=2)) )
          !--- First assign the interior points (excluding overlap lons)
          ocn_mask(2:nlon_o+1,:) = tmp_mask
          !--- Then assign the extra longitudes
          ocn_mask(1,:) = tmp_mask(nlon_o,:)
          ocn_mask(nlon_o+2,:) = tmp_mask(1,:)
        else
          write(6,*)"get_ocn_mask: Shape mismatch assigning ocn_mask with overlap lon"
          write(6,*)"    expected shape = ",nlon_o,nlat_o
          write(6,*)"       found shape = ",shape(tmp_mask)
          call xit("get_ocn_mask",-2)
        endif
      else
        allocate( ocn_mask(size(tmp_mask,dim=1),size(tmp_mask,dim=2)) )
        ocn_mask = nint(tmp_mask)
      endif
      deallocate(tmp_mask)

    end subroutine get_ocn_mask

    !***************************************************************************
    !--- Extract the 2D area values from the ocn grid definition
    !***************************************************************************
    subroutine get_ocn_area(ocn_area, overlap)
      real(kind=8), pointer, intent(inout) :: ocn_area(:,:)
      logical, intent(in), optional :: overlap

      !--- Local
      logical :: add_overlap
      real(kind=8), pointer :: tmp_area(:,:) => null()

      if ( present(overlap) ) then
        add_overlap = overlap
      else
        add_overlap = .false.
      endif

      if (id_ocn_grid < 0 ) then
        write(6,*)"get_ocn_area: OCN grid is not defined"
        call xit("get_ocn_area",-1)
      endif

      !--- Retrieve OCN cell center area values
      call get_grid_item(tmp_area, id_ocn_grid, "center", "area")

      if ( associated(ocn_area) ) deallocate(ocn_area)
      if (add_overlap) then
        if ( nlon_o == size(tmp_area,dim=1) ) then
          !--- Add the overlap longitude
          allocate( ocn_area(nlon_o+2,size(tmp_area,dim=2)) )
          !--- First assign the interior points (excluding overlap lons)
          ocn_area(2:nlon_o+1,:) = tmp_area
          !--- Then assign the extra longitudes
          ocn_area(1,:) = tmp_area(nlon_o,:)
          ocn_area(nlon_o+2,:) = tmp_area(1,:)
        else
          write(6,*)"get_ocn_area: Shape mismatch assigning ocn_area with overlap lon"
          write(6,*)"    expected shape = ",nlon_o,nlat_o
          write(6,*)"       found shape = ",shape(tmp_area)
          call xit("get_ocn_area",-2)
        endif
      else
        allocate( ocn_area(size(tmp_area,dim=1),size(tmp_area,dim=2)) )
        ocn_area = nint(tmp_area)
      endif
      deallocate(tmp_area)

    end subroutine get_ocn_area

    !***************************************************************************
    !--- Extract 2D integer*4 coordinate values from a defined grid
    !***************************************************************************
    subroutine get_grid_coord_2di4(coord, gridID, stagger, cdim)
      !--- The coordinate values are returned in coord
      !--- coord will be destroyed and recreated below
      integer(kind=4), pointer, intent(inout) :: coord(:,:)

      !--- The internal ID of the grid for which coords are required
      integer, intent(in) :: gridID

      !--- The stagger location (center or corner)
      character(*), intent(in) :: stagger

      !--- The grid dimension for which coords are required (1 or 2)
      integer, intent(in) :: cdim

      !--- Local
      real(kind=8), pointer :: coord_2dr8(:,:) => null()

      !--- Read coordinate data into a real*8 array
      call get_grid_coord_2dr8(coord_2dr8, gridID, stagger, cdim)

      !--- Assign this real*8 data to the input array
      if ( associated(coord) ) deallocate(coord)
      allocate( coord(size(coord_2dr8,dim=1),size(coord_2dr8,dim=2)) )
      coord = coord_2dr8

      !--- Clean up
      deallocate(coord_2dr8)

    end subroutine get_grid_coord_2di4

    !***************************************************************************
    !--- Extract 2D integer*4 coordinate values from a defined grid
    !***************************************************************************
    subroutine get_grid_coord_2di8(coord, gridID, stagger, cdim)
      !--- The coordinate values are returned in coord
      !--- coord will be destroyed and recreated below
      integer(kind=8), pointer, intent(inout) :: coord(:,:)

      !--- The internal ID of the grid for which coords are required
      integer, intent(in) :: gridID

      !--- The stagger location (center or corner)
      character(*), intent(in) :: stagger

      !--- The grid dimension for which coords are required (1 or 2)
      integer, intent(in) :: cdim

      !--- Local
      real(kind=8), pointer :: coord_2dr8(:,:) => null()

      !--- Read coordinate data into a real*8 array
      call get_grid_coord_2dr8(coord_2dr8, gridID, stagger, cdim)

      !--- Assign this real*8 data to the input array
      if ( associated(coord) ) deallocate(coord)
      allocate( coord(size(coord_2dr8,dim=1),size(coord_2dr8,dim=2)) )
      coord = coord_2dr8

      !--- Clean up
      deallocate(coord_2dr8)

    end subroutine get_grid_coord_2di8

    !***************************************************************************
    !--- Extract 2D real*4 coordinate values from a defined grid
    !***************************************************************************
    subroutine get_grid_coord_2dr4(coord, gridID, stagger, cdim)
      !--- The coordinate values are returned in coord
      !--- coord will be destroyed and recreated below
      real(kind=4), pointer, intent(inout) :: coord(:,:)

      !--- The internal ID of the grid for which coords are required
      integer, intent(in) :: gridID

      !--- The stagger location (center or corner)
      character(*), intent(in) :: stagger

      !--- The grid dimension for which coords are required (1 or 2)
      integer, intent(in) :: cdim

      !--- Local
      real(kind=8), pointer :: coord_2dr8(:,:) => null()

      !--- Read coordinate data into a real*8 array
      call get_grid_coord_2dr8(coord_2dr8, gridID, stagger, cdim)

      !--- Assign this real*8 data to the input array
      if ( associated(coord) ) deallocate(coord)
      allocate( coord(size(coord_2dr8,dim=1),size(coord_2dr8,dim=2)) )
      coord = coord_2dr8

      !--- Clean up
      deallocate(coord_2dr8)

    end subroutine get_grid_coord_2dr4

    !***************************************************************************
    !--- Extract 2D real*8 coordinate values from a defined grid
    !***************************************************************************
    subroutine get_grid_coord_2dr8(coord, gridID, stagger, cdim)
      !--- The coordinate values are returned in coord
      !--- coord will be destroyed and recreated below
      real(kind=8), pointer, intent(inout) :: coord(:,:)

      !--- The internal ID of the grid for which coords are required
      integer, intent(in) :: gridID

      !--- The stagger location (center or corner)
      character(*), intent(in) :: stagger

      !--- The grid dimension for which coords are required (1 or 2)
      integer, intent(in) :: cdim

      !--- Local
      character(64) :: method = "get_grid_coord_2dr8"
      integer(ESMF_KIND_I) :: lnum, rc, coordDim
      real(ESMF_KIND_R8), pointer :: coordPtr(:,:) => null()
      type (ESMF_StaggerLoc) :: staggerLoc
      character(256) :: strng

      if ( gridID < 1 .or. gridID > n_grid_list ) then
        write(6,*)"get_grid_coord_2dr8: Grid ",gridID," is not defined"
        call xit("get_grid_coord_2dr8",-1)
      endif

      strng = " "
      strng = lowerc( trim(adjustl(stagger)) )
      select case ( trim(strng) )
        case ("center", "cent", "mid")
          staggerLoc = ESMF_STAGGERLOC_CENTER
        case ("corner", "corn")
          staggerLoc = ESMF_STAGGERLOC_CORNER
        case ("edge1")
          staggerLoc = ESMF_STAGGERLOC_EDGE1
        case ("edge2")
          staggerLoc = ESMF_STAGGERLOC_EDGE2
        case default
          write(6,*)"get_grid_coord_2dr8: Invalid stagger location ",trim(stagger)
          call xit("get_grid_coord_2dr8",-2)
      end select

      !--- Use a local variable of the correct type to pass to ESMF_GridGetCoord
      coordDim = cdim

      !--- Retrieve a pointer to the coordinate values
      nullify(coordPtr)
      call ESMF_GridGetCoord(grid_list(gridID), coordDim=coordDim, localDE=zero, &
                             staggerLoc=staggerLoc, farrayPtr=coordPtr, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to retrieve coord values", line=lnum, method=method)

      if ( .not. associated(coordPtr) ) then
        write(6,*)"get_grid_coord_2dr8: Unable to retrieve coordinate data"
        call xit("get_grid_coord_2dr8",-3)
      endif

      !xxx if ( associated(coord) ) deallocate(coord)
      allocate( coord(size(coordPtr,dim=1),size(coordPtr,dim=2)) )
      coord = coordPtr
      nullify(coordPtr)

    end subroutine get_grid_coord_2dr8

    !***************************************************************************
    !--- Extract mask or area values from a defined grid
    !***************************************************************************
    subroutine get_grid_item(item_data, gridID, stagger, item)
      !--- The item values are returned in item_data
      !--- item_data will be destroyed and recreated below
      real(kind=8), pointer, intent(out) :: item_data(:,:)

      !--- The internal ID of the grid associated with this item
      integer, intent(in) :: gridID

      !--- The stagger location (center or corner)
      character(*), intent(in) :: stagger

      !--- The grid item that is requested (mask or area)
      character(*), intent(in) :: item

      !--- Local
      character(64) :: method = "get_grid_item"
      integer(ESMF_KIND_I) :: lnum, rc, rank
      type(ESMF_TypeKind_Flag) :: typekind
      real(ESMF_KIND_R4), pointer    :: item_2dr4(:,:) => null()
      real(ESMF_KIND_R8), pointer    :: item_2dr8(:,:) => null()
      integer(ESMF_KIND_I4), pointer :: item_2di4(:,:) => null()
      integer(ESMF_KIND_I8), pointer :: item_2di8(:,:) => null()
      type (ESMF_Array) :: itemArray
      type (ESMF_StaggerLoc) :: staggerLoc
      character(256) :: strng

      if ( gridID < 1 .or. gridID > n_grid_list ) then
        write(6,*)"get_grid_item: Grid ",gridID," is not defined"
        call xit("get_grid_item",-1)
      endif

      !--- Identify the grid stagger location
      strng = " "
      strng = lowerc( trim(adjustl(stagger)) )
      select case ( trim(strng) )
        case ("center", "cent", "mid")
          staggerLoc = ESMF_STAGGERLOC_CENTER
        case ("corner", "corn")
          staggerLoc = ESMF_STAGGERLOC_CORNER
        case ("edge1")
          staggerLoc = ESMF_STAGGERLOC_EDGE1
        case ("edge2")
          staggerLoc = ESMF_STAGGERLOC_EDGE2
        case default
          write(6,*)"get_grid_item: Invalid stagger location ",trim(stagger)
          call xit("get_grid_item",-2)
      end select

      !--- Identify the item that is to be returned
      strng = " "
      strng = lowerc( trim(adjustl(item)) )
      select case ( trim(strng) )
        case ("mask")
          !--- Retrieve an ESMF array containing the grid mask
          call ESMF_GridGetItem(grid_list(gridID), &
                                staggerLoc=staggerLoc, &
                                itemflag=ESMF_GRIDITEM_MASK, &
                                array=itemArray, rc=rc)
          lnum = __LINE__ ; call handle_ESMF_err(rc, &
              "Cannot retrieve grid mask", line=lnum, method=method)

        case ("area")
          !--- Retrieve an ESMF array containing the grid area
          call ESMF_GridGetItem(grid_list(gridID), &
                                staggerLoc=staggerLoc, &
                                itemflag=ESMF_GRIDITEM_AREA, &
                                array=itemArray, rc=rc)
          lnum = __LINE__ ; call handle_ESMF_err(rc, &
              "Cannot retrieve grid area", line=lnum, method=method)

        case default
          write(6,*)"get_grid_item: Invalid item ",trim(item)
          call xit("get_grid_item",-2)

      end select

      !--- Determine the type, kind and rank of the item array
      call ESMF_ArrayGet(itemArray, rank=rank, typeKind=typeKind, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to retrieve rank of dst_grid coord array", line=lnum, method=method)

      if ( rank /= 2 ) then
        write(6,*)"get_grid_item: Expecting rank=2 but found rank=",rank," for ",trim(item)
        call xit("get_grid_item",-2)
      endif

      !--- Destroy any existing item data
      if ( associated(item_data) ) deallocate(item_data)

      if ( typeKind == ESMF_TYPEKIND_I4 ) then
        call ESMF_ArrayGet(itemArray, farrayPtr=item_2di4, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve pointer to "//trim(item), line=lnum, method=method)
        if ( .not. associated(item_2di4) ) then
          write(6,*)"get_grid_item: Unable to retreive 2DI4 item ",trim(item)
          call xit("get_grid_item",-2)
        endif
        allocate( item_data(size(item_2di4,dim=1),size(item_2di4,dim=2)) )
        item_data = item_2di4
        nullify(item_2di4)

      else if ( typeKind == ESMF_TYPEKIND_I8 ) then
        call ESMF_ArrayGet(itemArray, farrayPtr=item_2di8, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve pointer to "//trim(item), line=lnum, method=method)
        if ( .not. associated(item_2di8) ) then
          write(6,*)"get_grid_item: Unable to retreive 2DI8 item ",trim(item)
          call xit("get_grid_item",-2)
        endif
        allocate( item_data(size(item_2di8,dim=1),size(item_2di8,dim=2)) )
        item_data = item_2di8
        nullify(item_2di8)

      else if ( typeKind == ESMF_TYPEKIND_R4 ) then
        call ESMF_ArrayGet(itemArray, farrayPtr=item_2dr4, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve pointer to "//trim(item), line=lnum, method=method)
        if ( .not. associated(item_2dr4) ) then
          write(6,*)"get_grid_item: Unable to retreive 2DR4 item ",trim(item)
          call xit("get_grid_item",-2)
        endif
        allocate( item_data(size(item_2dr4,dim=1),size(item_2dr4,dim=2)) )
        item_data = item_2dr4
        nullify(item_2dr4)

      else if ( typeKind == ESMF_TYPEKIND_R8 ) then
        call ESMF_ArrayGet(itemArray, farrayPtr=item_2dr8, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve pointer to "//trim(item), line=lnum, method=method)
        if ( .not. associated(item_2dr8) ) then
          write(6,*)"get_grid_item: Unable to retreive 2DR8 item ",trim(item)
          call xit("get_grid_item",-2)
        endif
        allocate( item_data(size(item_2dr8,dim=1),size(item_2dr8,dim=2)) )
        item_data = item_2dr8
        nullify(item_2dr8)

      else
        write(6,*)"get_grid_item: Invalid typeKind"
        call xit("get_grid_item",-2)
      endif

    end subroutine get_grid_item

    !***************************************************************************
    !--- Create a netcd file containing AGCM grid information
    !--- This will contain the coordinate, fractional land mask and cell area
    !--- data that is available at run time
    !***************************************************************************
    subroutine create_atm_grid_file(fname)
      character(*), intent(in) :: fname

      !--- Local
      integer(kind=4) :: fid, vid(20), dimid(10)
      real(kind=8), pointer :: mid_lon(:,:) => null(), mid_lat(:,:) => null()
      real(kind=8), pointer :: corn_lon(:,:) => null(), corn_lat(:,:) => null()
      real(kind=8), pointer :: cell_area(:,:) => null(), fland(:,:) => null()
      real(kind=8), pointer :: wrk2d(:,:) => null()
      real(kind=8) :: avg_area
      integer :: nlon, nlat, nlon_corn, nlat_corn

      if ( len_trim(fname) < 1 ) then
        write(6,*)"create_atm_grid_file: No file name was supplied"
        call xit("create_atm_grid_file",-1)
      endif

      !--- Data from the internal grid descriptions found in grid_list

      !--- Retrieve ATM cell center lon (dim 1) coordinate values
      call get_grid_coord(mid_lon, id_atm_grid, "center", 1)

      !--- Retrieve ATM cell center lat (dim 2) coordinate values
      call get_grid_coord(mid_lat, id_atm_grid, "center", 2)

      !--- Retrieve ATM cell corner lon (dim 1) coordinate values
      call get_grid_coord(corn_lon, id_atm_grid, "corner", 1)

      !--- Retrieve ATM cell corner lat (dim 2) coordinate values
      call get_grid_coord(corn_lat, id_atm_grid, "corner", 2)

      !--- Retrieve ATM grid cell area values
      call get_grid_item(cell_area, id_atm_grid, "center", "area")
      avg_area = sum(cell_area)/real(size(cell_area),kind=8)

      nlon = size(mid_lon,dim=1)
      nlat = size(mid_lon,dim=2)
      nlon_corn = size(corn_lon,dim=1)
      nlat_corn = size(corn_lon,dim=2)

      !--- Open the file and define dimensions
      call ncdf_open(fid, fname, "create")
      call ncdf_add_dim(fid, dimid(1), "nlon", nlon)
      call ncdf_add_dim(fid, dimid(2), "nlat", nlat)
      call ncdf_add_dim(fid, dimid(3), "nlon_corn", nlon_corn)
      call ncdf_add_dim(fid, dimid(4), "nlat_corn", nlat_corn)
      call ncdf_add_dim(fid, dimid(5), "ccrecl", (nlon+1)*nlat)

      !--- Define the earth radius (according to the AGCM) as a global attribute
      call ncdf_add_att(fid, NF90_GLOBAL, "AGCM_earth_radius_in_meters", atm_earth_radius)

      !--- Define and assign center longitude
      call ncdf_add_var(fid, vid(1), "mid_lon", dimid(1:2), vals=mid_lon)
      call ncdf_add_att(fid, vid(1), "long_name", "AGCM Cell Center Longitude")
      call ncdf_add_att(fid, vid(1), "units", "degrees_east")

      !--- Define and assign center latitude
      call ncdf_add_var(fid, vid(2), "mid_lat", dimid(1:2), vals=mid_lat)
      call ncdf_add_att(fid, vid(2), "long_name", "AGCM Cell Center Latitude")
      call ncdf_add_att(fid, vid(2), "units", "degrees_north")

      !--- Define and assign corner longitude
      call ncdf_add_var(fid, vid(3), "corn_lon", dimid(3:4), vals=corn_lon)
      call ncdf_add_att(fid, vid(3), "long_name", "AGCM Cell Corner Longitude")
      call ncdf_add_att(fid, vid(3), "units", "degrees_east")

      !--- Define and assign corner latitude
      call ncdf_add_var(fid, vid(4), "corn_lat", dimid(3:4), vals=corn_lat)
      call ncdf_add_att(fid, vid(4), "long_name", "AGCM Cell Corner Latitude")
      call ncdf_add_att(fid, vid(4), "units", "degrees_north")

      !--- Define and assign cell area values
      call ncdf_add_var(fid, vid(5), "area_weights", dimid(1:2), vals=cell_area)
      call ncdf_add_att(fid, vid(5), "long_name", "AGCM Grid Cell Area Weights")
      call ncdf_add_att(fid, vid(5), "units", "1")
      call ncdf_add_att(fid, vid(5), "comment", "Cell areas used to calculate conservative interpolation weights")

      if ( associated(atm_fland) )then
        !--- Define and assign fractional land mask
        if ( size(atm_fland) /= (nlon+1)*nlat ) then
          write(6,*)"create_atm_grid_file: Unexpected size of atm_fland"
          write(6,*)"    expected ",(nlon+1)*nlat
          write(6,*)"      found ",size(atm_fland)
          call xit("create_atm_grid_file",-1)
        endif

        !--- Reformat atm_fland as a 2D array with the overlap longitude removed
        if ( associated(wrk2d) ) deallocate(wrk2d)
        allocate( wrk2d(nlon+1,nlat) )
        wrk2d = reshape(atm_fland, (/nlon+1,nlat/))
        if ( associated(fland) ) deallocate(fland)
        allocate( fland(nlon,nlat) )
        fland(:,:) = wrk2d(1:nlon,:)
        deallocate(wrk2d)
        call ncdf_add_var(fid, vid(7), "atm_fland_2d", dimid(1:2), vals=fland)
        call ncdf_add_att(fid, vid(7), "long_name", "AGCM Fractional Land Mask (0=ocean, 1=land)")
        call ncdf_add_att(fid, vid(7), "units", "1")
      endif

      !--- Data sent from the agcm to the coupler at run time

      if ( associated(atm_grid_cell_area) )then
        !--- Define and assign cell area values
        call ncdf_add_var(fid, vid(7), "atm_grid_cell_area", dimid(1:2), vals=atm_grid_cell_area)
        call ncdf_add_att(fid, vid(7), "long_name", "AGCM Grid Cell Area")
        call ncdf_add_att(fid, vid(7), "units", "m2")
        call ncdf_add_att(fid, vid(7), &
             "comment", "Dimensional values sent from the AGCM to the coupler at run time")
      endif

      if ( associated(atm_fland) )then
        !--- Write the 1D atm_fland array that was received from the AGCM "as is"
        call ncdf_add_var(fid, vid(7), "atm_fland", (/dimid(5)/), vals=atm_fland)
        call ncdf_add_att(fid, vid(7), "long_name", "AGCM Fractional Land Mask (0=ocean, 1=land)")
        call ncdf_add_att(fid, vid(7), "units", "1")
        call ncdf_add_att(fid, vid(7), &
             "comment", "Sent from the AGCM to the coupler at run time (contains overlap lon)")
      endif

      if ( associated(atm_lon) )then
        call ncdf_add_var(fid, vid(7), "atm_lon", (/dimid(1)/))
        call ncdf_add_att(fid, vid(7), "long_name", "AGCM Longitude")
        call ncdf_add_att(fid, vid(7), "units", "degrees_east")
        call ncdf_add_att(fid, vid(7), "comment", "Sent from the AGCM to the coupler at run time")
        call ncdf_write_var(fid, vid(7), vals=atm_lon)
      endif

      if ( associated(atm_lat) )then
        call ncdf_add_var(fid, vid(7), "atm_lat", (/dimid(2)/))
        call ncdf_add_att(fid, vid(7), "long_name", "AGCM Latitude")
        call ncdf_add_att(fid, vid(7), "units", "degrees_north")
        call ncdf_add_att(fid, vid(7), "comment", "Sent from the AGCM to the coupler at run time")
        call ncdf_write_var(fid, vid(7), vals=atm_lat)
      endif

      if ( associated(atm_gwt) )then
        call ncdf_add_var(fid, vid(7), "atm_gwt", (/dimid(2)/))
        call ncdf_add_att(fid, vid(7), "long_name", "Gaussian Weights")
        call ncdf_add_att(fid, vid(7), "units", "1")
        call ncdf_add_att(fid, vid(7), "comment", "Sent from the AGCM to the coupler at run time")
        call ncdf_write_var(fid, vid(7), vals=atm_gwt)
      endif

      !--- Close the file
      call ncdf_close(fid)

      if ( associated(mid_lon) )   deallocate(mid_lon)
      if ( associated(mid_lat) )   deallocate(mid_lat)
      if ( associated(corn_lon) )  deallocate(corn_lon)
      if ( associated(corn_lat) )  deallocate(corn_lat)
      if ( associated(cell_area) ) deallocate(cell_area)
      if ( associated(fland) )     deallocate(fland)

    end subroutine create_atm_grid_file

    !***************************************************************************
    !--- Create a netcd file containing NEMO grid information
    !--- This will contain the coordinate, mask and cell area data that is
    !--- available at run time
    !***************************************************************************
    subroutine create_ocn_grid_file(fname)
      character(*), intent(in) :: fname

      !--- Local
      integer(kind=4) :: fid, vid(20), dimid(10)
      real(kind=8), pointer :: mid_lon(:,:) => null(), mid_lat(:,:) => null()
      real(kind=8), pointer :: corn_lon(:,:) => null(), corn_lat(:,:) => null()
      real(kind=8), pointer :: cell_area(:,:) => null(), wrk2d(:,:) => null()
      integer, pointer :: cell_mask(:,:) => null()
      real(kind=8) :: avg_area
      integer :: nlon, nlat, nlon_corn, nlat_corn

      if ( len_trim(fname) < 1 ) then
        write(6,*)"create_ocn_grid_file: No file name was supplied"
        call xit("create_ocn_grid_file",-1)
      endif

      !--- Retrieve OCN cell center lon (dim 1) coordinate values
      call get_grid_coord(mid_lon, id_ocn_grid, "center", 1)

      !--- Retrieve OCN cell center lat (dim 2) coordinate values
      call get_grid_coord(mid_lat, id_ocn_grid, "center", 2)

      !--- Retrieve OCN cell corner lon (dim 1) coordinate values
      call get_grid_coord(corn_lon, id_ocn_grid, "corner", 1)

      !--- Retrieve OCN cell corner lat (dim 2) coordinate values
      call get_grid_coord(corn_lat, id_ocn_grid, "corner", 2)

      !--- Retrieve OCN grid cell area values
      call get_grid_item(cell_area, id_ocn_grid, "center", "area")
      avg_area = sum(cell_area)/real(size(cell_area),kind=8)

      call get_ocn_mask(cell_mask)

      nlon = size(mid_lon,dim=1)
      nlat = size(mid_lon,dim=2)
      nlon_corn = size(corn_lon,dim=1)
      nlat_corn = size(corn_lon,dim=2)

      !--- Open the file and define dimensions
      call ncdf_open(fid, fname, "create")
      call ncdf_add_dim(fid, dimid(1), "nlon", nlon)
      call ncdf_add_dim(fid, dimid(2), "nlat", nlat)
      call ncdf_add_dim(fid, dimid(3), "nlon_corn", nlon_corn)
      call ncdf_add_dim(fid, dimid(4), "nlat_corn", nlat_corn)
      call ncdf_add_dim(fid, dimid(5), "jpiglo", nemo_jpiglo)
      call ncdf_add_dim(fid, dimid(6), "jpjglo", nemo_jpjglo)

      !--- Define the earth radius (according to the AGCM) as a global attribute
      call ncdf_add_att(fid, NF90_GLOBAL, "AGCM_earth_radius_in_meters", atm_earth_radius)

      !--- Define and assign center longitude
      call ncdf_add_var(fid, vid(1), "mid_lon", dimid(1:2), vals=mid_lon)
      call ncdf_add_att(fid, vid(1), "long_name", "NEMO T Grid Cell Center Longitude")
      call ncdf_add_att(fid, vid(1), "units", "degrees_east")

      !--- Define and assign center latitude
      call ncdf_add_var(fid, vid(2), "mid_lat", dimid(1:2), vals=mid_lat)
      call ncdf_add_att(fid, vid(2), "long_name", "NEMO T Grid Cell Center Latitude")
      call ncdf_add_att(fid, vid(2), "units", "degrees_north")

      !--- Define and assign corner longitude
      call ncdf_add_var(fid, vid(3), "corn_lon", dimid(3:4), vals=corn_lon)
      call ncdf_add_att(fid, vid(3), "long_name", "NEMO Grid Cell Corner Longitude")
      call ncdf_add_att(fid, vid(3), "units", "degrees_east")

      !--- Define and assign corner latitude
      call ncdf_add_var(fid, vid(4), "corn_lat", dimid(3:4), vals=corn_lat)
      call ncdf_add_att(fid, vid(4), "long_name", "NEMO Grid Cell Corner Latitude")
      call ncdf_add_att(fid, vid(4), "units", "degrees_north")

      !--- Define and assign cell mask values (1=ocean, 0=land)
      if ( associated(wrk2d) ) deallocate(wrk2d)
      allocate( wrk2d(nlon,nlat) )
      wrk2d = cell_mask
      call ncdf_add_var(fid, vid(5), "mask", dimid(1:2), vals=wrk2d)
      call ncdf_add_att(fid, vid(5), "long_name", "NEMO T Grid Cell Mask (1=ocean, 0=land)")
      call ncdf_add_att(fid, vid(5), "units", "1")

      !--- Define and assign cell area weights values
      call ncdf_add_var(fid, vid(6), "area_weights", dimid(1:2), vals=cell_area)
      call ncdf_add_att(fid, vid(6), "long_name", "NEMO Grid Cell Area Weights")
      call ncdf_add_att(fid, vid(6), "units", "1")
      call ncdf_add_att(fid, vid(6), "comment", "Cell areas used to calculate conservative interpolation weights")

      if ( associated(nemo_glamt) ) then
        call ncdf_add_var(fid, vid(7), "glamt", dimid(5:6), vals=nemo_glamt)
        call ncdf_add_att(fid, vid(7), "comment", "Sent from NEMO to the coupler at run time")
      endif

      if ( associated(nemo_gphit) ) then
        call ncdf_add_var(fid, vid(7), "gphit", dimid(5:6), vals=nemo_gphit)
        call ncdf_add_att(fid, vid(7), "comment", "Sent from NEMO to the coupler at run time")
      endif

      if ( associated(nemo_glamf) ) then
        call ncdf_add_var(fid, vid(7), "glamf", dimid(5:6), vals=nemo_glamf)
        call ncdf_add_att(fid, vid(7), "comment", "Sent from NEMO to the coupler at run time")
      endif

      if ( associated(nemo_gphif) ) then
        call ncdf_add_var(fid, vid(7), "gphif", dimid(5:6), vals=nemo_gphif)
        call ncdf_add_att(fid, vid(7), "comment", "Sent from NEMO to the coupler at run time")
      endif

      if ( associated(nemo_tmask) ) then
        call ncdf_add_var(fid, vid(7), "tmask", dimid(5:6), vals=nemo_tmask)
        call ncdf_add_att(fid, vid(7), "comment", "Sent from NEMO to the coupler at run time")
      endif

      if ( associated(nemo_e1t) ) then
        call ncdf_add_var(fid, vid(7), "e1t", dimid(5:6), vals=nemo_e1t)
        call ncdf_add_att(fid, vid(7), "comment", "Sent from NEMO to the coupler at run time")
      endif

      if ( associated(nemo_e2t) ) then
        call ncdf_add_var(fid, vid(7), "e2t", dimid(5:6), vals=nemo_e2t)
        call ncdf_add_att(fid, vid(7), "comment", "Sent from NEMO to the coupler at run time")
      endif

      !--- Close the file
      call ncdf_close(fid)

      if ( associated(mid_lon) )   deallocate(mid_lon)
      if ( associated(mid_lat) )   deallocate(mid_lat)
      if ( associated(corn_lon) )  deallocate(corn_lon)
      if ( associated(corn_lat) )  deallocate(corn_lat)
      if ( associated(cell_area) ) deallocate(cell_area)
      if ( associated(cell_mask) ) deallocate(cell_mask)
      if ( associated(wrk2d) )     deallocate(wrk2d)

    end subroutine create_ocn_grid_file

    !***************************************************************************
    !--- Read coordinate information from a netcd file containing AGCM grid
    !--- data that was created by create_atm_grid_file
    !***************************************************************************
    subroutine read_atm_grid_file(fname)
      character(*), intent(in) :: fname

      !--- Local
      integer(kind=4) :: fid, ncstat

      if ( len_trim(fname) < 1 ) then
        write(6,*)"read_atm_grid_file: No file name was supplied"
        call xit("read_atm_grid_file",-1)
      endif

      !--- Open the file for reading
      call ncdf_open(fid, fname, "read")

      !--- Read selected grid information from the file
      !--- The following calls to ncdf_read_var will destroy
      !--- and then recreate each data array
      call ncdf_read_var(fid, "atm_lon",   atm_lon)
      call ncdf_read_var(fid, "atm_lat",   atm_lat)
      call ncdf_read_var(fid, "atm_gwt",   atm_gwt)
      call ncdf_read_var(fid, "atm_fland", atm_fland)
      call ncdf_read_var(fid, "atm_grid_cell_area", atm_grid_cell_area)

      ncstat = nf90_get_att(fid,NF90_GLOBAL,'AGCM_earth_radius_in_meters',atm_earth_radius)
      call nc_error_handler(fid, ncstat, 1, &
          msg="read_atm_grid_file: Reading global att AGCM_earth_radius_in_meters")

      !--- Close the file
      call ncdf_close(fid)

      !--- Define AGCM dimension sizes from the shape of the lon and lat arrays
      atm_nlon = size(atm_lon)
      atm_nlat = size(atm_lat)

      write(6,*)"read_atm_grid_file: Read AGCM grid information from ",trim(fname), &
          "  grid shape = ",atm_nlon,atm_nlat

    end subroutine read_atm_grid_file

    !***************************************************************************
    !--- Read coordinate information from a netcd file containing NEMO grid
    !--- data that was created by create_ocn_grid_file
    !***************************************************************************
    subroutine read_ocn_grid_file(fname)
      character(*), intent(in) :: fname

      !--- Local
      integer(kind=4) :: fid, ncstat

      if ( len_trim(fname) < 1 ) then
        write(6,*)"read_ocn_grid_file: No file name was supplied"
        call xit("read_ocn_grid_file",-1)
      endif

      !--- Open the file for reading
      call ncdf_open(fid, fname, "read")

      !--- Read selected grid information from the file
      !--- The following calls to ncdf_read_var will destroy
      !--- and then recreate the nemo_* grid data arrays
      call ncdf_read_var(fid, "glamt", nemo_glamt)
      call ncdf_read_var(fid, "gphit", nemo_gphit)
      call ncdf_read_var(fid, "glamf", nemo_glamf)
      call ncdf_read_var(fid, "gphif", nemo_gphif)
      call ncdf_read_var(fid, "e1t",   nemo_e1t)
      call ncdf_read_var(fid, "e2t",   nemo_e2t)
      call ncdf_read_var(fid, "tmask", nemo_tmask)

      if ( atm_earth_radius < 0.0 ) then
        !--- Only read the earth radius if it is not already defined
        ncstat = nf90_get_att(fid,NF90_GLOBAL,'AGCM_earth_radius_in_meters',atm_earth_radius)
        call nc_error_handler(fid, ncstat, 1, &
            msg="read_ocn_grid_file: Reading global att AGCM_earth_radius_in_meters")
      endif

      !--- Close the file
      call ncdf_close(fid)

      !--- Define NEMO global dimension sizes from the shape of these arrays
      nemo_jpiglo = size(nemo_glamt,dim=1)
      nemo_jpjglo = size(nemo_glamt,dim=2)

      write(6,*)"read_ocn_grid_file: Read NEMO grid information from ",trim(fname), &
          "  grid shape = ",nemo_jpiglo,nemo_jpjglo

    end subroutine read_ocn_grid_file

    !***************************************************************************
    !--- Copy one cpl_grid_t variable to another cpl_grid_t variable
    !***************************************************************************
    subroutine copy_cpl_grid_t(in, out)
      type(cpl_grid_t), intent( in) :: in
      type(cpl_grid_t), intent(out) :: out
      out % title     = in % title
      out % file_name = in % file_name
      out % nx        = in % nx
      out % ny        = in % ny
      out % size      = in % size
      out % edges     = in % edges
      if ( associated( in%mask ) ) then
        if ( associated(out%mask) ) deallocate(out%mask)
        allocate(out%mask(size(in%mask,dim=1),size(in%mask,dim=2)))
        out % mask = in % mask
      else
        out % mask => null()
      endif
      if ( associated( in%grid ) ) then
        !--- out%grid is an alias for in%grid
        out % grid = in % grid
      else
        out % grid => null()
      endif
    end subroutine copy_cpl_grid_t

    !***************************************************************************
    !--- Error handling for ESMF routines
    !***************************************************************************
    subroutine handle_ESMF_err(rc, txt, line, method)

      integer(ESMF_KIND_I), intent(in) :: rc
      character(*),         intent(in) :: txt
      integer(ESMF_KIND_I), intent(in), optional :: line
      character(*),         intent(in), optional :: method

      !--- Local
      logical :: found

      if ( present(line) .and. present(method) ) then
        found = ESMF_LogFoundError(rc, msg=trim(txt), line=line, method=trim(method))
      else if ( present(line) ) then
        found = ESMF_LogFoundError(rc, msg=trim(txt), line=line)
      else if ( present(method) ) then
        found = ESMF_LogFoundError(rc, msg=trim(txt), method=trim(method))
      else
        found = ESMF_LogFoundError(rc, msg=trim(txt))
      endif
      if ( found ) then
        if ( len(trim(txt)) > 0 ) then
          if ( present(method) .and. present(line) ) then
            write(6,'(3a,i8,2x,a)') " ## EE ## ",trim(method),":",line,trim(txt)
          else if ( present(method) ) then
            write(6,'(2a,2x,a)') " ## EE ## ",trim(method),trim(txt)
          else if ( present(line) ) then
            write(6,'(a,i8,a)') " ## EE ## ",line,trim(txt)
          else
            write(6,'(2a)') " ## EE ## ",trim(txt)
          endif
          call flush(6)
        endif
        write(6,*)'See the log file PET0.',trim(adjustl(default_log_file))
        call ESMF_LogFlush()
        call xit("ESMF_err",-1)
      endif

    end subroutine handle_ESMF_err

    !***************************************************************************
    !--- Return a string containing time on a given clock
    !***************************************************************************
    function get_time_string(this_clock, time_type) result(time_string)

      !--- The clock to read the current time from
      type(ESMF_Clock), intent(inout) :: this_clock

      !--- The type of time to return (default is current time)
      character(*), intent(in), optional :: time_type

      !--- The return value, time formatted in a string
      character(len=32) :: time_string

      type(ESMF_Time)         :: ref_time, curr_time
      type(ESMF_TimeInterval) :: sim_time
      integer(ESMF_KIND_I8) :: advance_count
      integer(ESMF_KIND_I4) :: year, month, day, hour, minute, second
      integer(ESMF_KIND_I)  :: rc
      character(32) :: which_time

      if ( present(time_type) ) then
        which_time = trim(time_type)
      else
        which_time = "current"
      endif

      !--- Get info about this clock
      call ESMF_ClockGet(this_clock, refTime=ref_time, currTime=curr_time, &
                         currSimTime=sim_time, advanceCount=advance_count, rc=rc)
      call handle_ESMF_err(rc, "get_time_string: Problem reading time.")

      select case ( trim(adjustl(which_time)) )
        case ("ref","reference","ref_time")
          !--- Return the reference time
          call ESMF_TimeGet(ref_time, yy=year, mm=month, dd=day, &
                                      h=hour, m=minute, s=second, rc=rc)
          call handle_ESMF_err(rc, "get_time_string: Problem getting reference time.")
        case ("curr","current","curr_time")
          !--- Return the current time
          call ESMF_TimeGet(curr_time, yy=year, mm=month, dd=day, &
                                       h=hour, m=minute, s=second, rc=rc)
          call handle_ESMF_err(rc, "get_time_string: Problem getting current time.")
        case ("sim","simulation","sim_time")
          !--- Return simulation time (currTime - refTime)
          call ESMF_TimeIntervalGet(sim_time, ref_time, yy=year, mm=month, d=day, &
                                      h=hour, m=minute, s=second, rc=rc)
          call handle_ESMF_err(rc, "get_time_string: Problem getting simulation time.")
        case default
          write(6,*)"get_time_string: Invalid time type ",trim(time_type)
          call flush(6)
          call ESMF_LogFlush()
          call ESMF_Finalize(endflag=ESMF_END_ABORT)
      end select

      !--- Format the output string
      time_string = " "
      write(time_string,'(i4.4,"-",i2.2,"-",i2.2,1x,i2.2,":",i2.2,":",i2.2)') &
            year,month,day,hour,minute,second

    end function get_time_string

    !***************************************************************************
    !--- Initialize ESMF
    !***************************************************************************
    subroutine init_esmf(mpi_group, logname)

      !--- The mpi communicator assigned to the group of the invoking routine
      integer, intent(in) :: mpi_group

      !--- The name of the default log file to be created
      character(*), intent(in) :: logname

      integer(ESMF_KIND_I) :: mpicomm, UnitLBound_esmf, UnitUBound_esmf
      integer(ESMF_KIND_I) :: rc
      integer(kind=impi) :: ierr
      logical(kind=impi) :: mpi_init_was_called

#ifdef use_mpi
      !--- mpi_init must be called prior to calling this routine
      call mpi_initialized( mpi_init_was_called, ierr )
      if ( .not. mpi_init_was_called ) then
        write(6,*)"init_esmf: MPI has not been initialized."
        call xit("init_esmf",-1)        
      endif

      !--- Initialize ESMF
      !---
      !--- defaultCalKind     ..The default calander to use in all ESMF routines
      !--- defaultlogfilename ..The name of the default log file used by ESMF routines
      !--- mpiCommunicator    ..Communicator defining the MPI group in which ESMF will run
      !--- ioUnitLBound       ..Lower bound for fortran unit numbers used by ESMF
      !--- ioUnitUBound       ..Upper bound for fortran unit numbers used by ESMF
      !--- logkindflag        ..Flags the use of none, one or many (1 per mpi task) log files

      !--- These parameters must be of default integer kind for the esmf lib
      mpicomm = mpi_group
      UnitLBound_esmf = 300
      UnitUBound_esmf = 349

      !--- Save the log file name for use by handle_ESMF_err and others
      default_log_file = trim(logname)

      if ( len(trim(logname)) > 0 ) then
        !--- Use the user supplied log file name
        call ESMF_Initialize ( defaultCalKind=ESMF_CALKIND_NOLEAP, &
                           defaultlogfilename=trim(logname), &
                           mpiCommunicator=mpicomm, &
                           ioUnitLBound=UnitLBound_esmf, &
                           ioUnitUBound=UnitUBound_esmf, &
                           logkindflag=ESMF_LOGKIND_MULTI, &
                           rc=rc)
      else
        !--- Use the ESMF default log file name
        call ESMF_Initialize ( defaultCalKind=ESMF_CALKIND_NOLEAP, &
                           mpiCommunicator=mpicomm, &
                           ioUnitLBound=UnitLBound_esmf, &
                           ioUnitUBound=UnitUBound_esmf, &
                           logkindflag=ESMF_LOGKIND_MULTI, &
                           rc=rc)
      endif
      call handle_ESMF_err(rc, "Problem with ESMF_Initialize")
#else
      write(6,*)"init_esmf: *** WARNING *** MPI is not being used."
#endif

    end subroutine init_esmf

    !***************************************************************************
    !--- Allocate and assign the global event array
    !--- This will define the set of events the coupler will handle
    !***************************************************************************
    subroutine define_events()

      integer :: idx1, idx2, nvars
      character(64) :: get_var, put_var
      character(32), allocatable :: tmp_uniq_names(:)
      integer, allocatable :: gids(:)
      character(1024) :: strng
      logical :: found
      type(cpl_vinfo_t) :: get_vinfo, put_vinfo

      !--- The global event array must be allocated and partially defined prior
      !--- to calling this routine. The following elements must be set
      !---    get_model ...the model from with the transfer will originate (atm, ocn, ice)
      !---    put_model ...the model in which the tranfer will be received (atm, ocn, ice)
      !---    get_var   ...the variable name associated with get_model
      !---    put_var   ...the variable name associated with put_model
      !--- get_var and put_var are names known to the coupler, meaning that they are
      !--- found in the internal array cpl_var_list (defined in com_cpl)
      !--- get_var and put_var may be the same, in which case no interpolation
      !--- is done in the coupler and the field is passed from get_model to
      !--- put_model without any modification.
      if ( nevents <= 0 ) then
        write(6,*)"define_events: The events array is not defined."
        call xit("define_events",-1)
      endif

      do idx1=1,nevents

        if ( len_trim(event(idx1)%get_model) < 1 ) then
          write(6,*)"define_events: get_model is not defined for event ",idx1
          call xit("define_events",-2)
        endif

        if ( len_trim(event(idx1)%put_model) < 1 ) then
          write(6,*)"define_events: put_model is not defined for event ",idx1
          call xit("define_events",-3)
        endif

        if ( len_trim(event(idx1)%get_var) < 1 ) then
          write(6,*)"define_events: get_var is not defined for event ",idx1
          call xit("define_events",-4)
        endif

        if ( len_trim(event(idx1)%put_var) < 1 ) then
          write(6,*)"define_events: put_var is not defined for event ",idx1
          call xit("define_events",-5)
        endif

        get_vinfo = find_cpl_vinfo( name=trim(adjustl(event(idx1)%get_var)) )
        put_vinfo = find_cpl_vinfo( name=trim(adjustl(event(idx1)%put_var)) )

        !--- Set a default regrid method if not set by the user
        if ( len_trim(event(idx1)%regrid_method) < 1 ) then
          event(idx1)%regrid_method = trim(adjustl( put_vinfo%regrid%method ))
        endif

        !--- Set a default regrid norm type if not set by the user
        if ( len_trim(event(idx1)%regrid_norm_opt) < 1 ) then
          event(idx1)%regrid_norm_opt = "fracarea"
        endif

        !--- Initialize the regrid_id to zero
        event(idx1)%regrid_id = 0

        !--- Assign get_master and put_master to be consistent with
        !--- the user supplied values for get_model and put_model
        strng = " "
        strng = lowerc( trim(adjustl(event(idx1)%get_model)) )
        select case ( trim(strng) )
          case ("cpl")
             event(idx1)%get_master = cpl_master
          case ("atm")
             event(idx1)%get_master = atm_master
          case ("ocn")
             event(idx1)%get_master = ocn_master
          case ("ice")
             event(idx1)%get_master = ice_master
          case default
            write(6,'(2a)')'define_events: Unknown get model ',trim(strng)
            call xit("define_events",-6)
        end select

        strng = " "
        strng = lowerc( trim(adjustl(event(idx1)%put_model)) )
        select case ( trim(strng) )
          case ("cpl")
             event(idx1)%put_master = cpl_master
          case ("atm")
             event(idx1)%put_master = atm_master
          case ("ocn")
             event(idx1)%put_master = ocn_master
          case ("ice")
             event(idx1)%put_master = ice_master
          case default
            write(6,'(2a)')'define_events: Unknown put model ',trim(strng)
            call xit("define_events",-6)
        end select

        !--- Never regrid a message transfer
        if ( event(idx1)%is_msg ) event(idx1)%regrid_id = 0

        !--- Never regrid a field that has the same name for get and put operations
        if ( trim(adjustl(event(idx1)%get_var)) .eq. &
             trim(adjustl(event(idx1)%put_var)) ) then
          event(idx1)%regrid_id = 0
        endif

        if (verbose > 0) then
          write(6,'(a,i3,2a)')' define_events: event ',idx1, &
                              '  action: ',trim(event(idx1)%action)
        endif

      enddo

      !--- Get the list of unique variable names for all transfer fields
      !--- The uniq_names array is defined in the module com_cpl but is
      !--- allocated and assigned here
      allocate( tmp_uniq_names(2*nevents) )
      nvars=0
      do idx1=1,nevents
        get_var = trim(adjustl(event(idx1)%get_var))
        put_var = trim(adjustl(event(idx1)%put_var))
        if (nvars == 0) then
          if ( .not. is_special_op(get_var) ) then
            nvars=1
            tmp_uniq_names(nvars) = trim(get_var)
          endif
          if ( .not. is_special_op(put_var) ) then
            if ( trim(get_var) .ne. trim(put_var) ) then
              nvars=nvars+1
              tmp_uniq_names(nvars) = trim(put_var)
            endif
          endif
        else
          !--- If get_var is not already in the list then append it
          found=.false.
          do idx2=1,nvars
            if ( trim(get_var) .eq. trim(tmp_uniq_names(idx2)) ) then
              !--- get_var is already in the list
              found=.true.
              exit
            endif
          enddo
          if ( .not. found .and. .not. is_special_op(get_var) ) then
            !--- Add get_var to the list
            nvars=nvars+1
            tmp_uniq_names(nvars) = trim(get_var)
          endif

          !--- If put_var is not already in the list then append it
          found=.false.
          do idx2=1,nvars
            if ( trim(put_var) .eq. trim(tmp_uniq_names(idx2)) ) then
              !--- put_var is already in the list
              found=.true.
              exit
            endif
          enddo
          if ( .not. found .and. .not. is_special_op(put_var) ) then
            !--- Add put_var to the list
            nvars=nvars+1
            tmp_uniq_names(nvars) = trim(put_var)
          endif
        endif
      enddo

      !---- Allocate and assign uniq_names
      if ( associated(uniq_names) ) deallocate(uniq_names)
      allocate( uniq_names(nvars), gids(nvars) )
      do idx1=1,nvars
        get_vinfo = find_cpl_vinfo( name=trim(adjustl(tmp_uniq_names(idx1))) )
        uniq_names(idx1)%name = trim(tmp_uniq_names(idx1))
        uniq_names(idx1)%grid = get_vinfo%grid
        !--- Set the group ID according to the grid type for this initialization
        !--- Additional name,gid pairs will need to be added later
        select case( trim(adjustl(get_vinfo%grid)) )
          case ("ocn")
            gids(idx1) = ocn_gid
          case ("atm")
            gids(idx1) = atm_gid
          case ("ice")
            gids(idx1) = ice_gid
          case ("msg")
            gids(idx1) = msg_gid
          case ("canom")
            gids(idx1) = canom_gid
          case default
            gids(idx1) = other_gid
        end select
      enddo

      deallocate( tmp_uniq_names, gids )

    end subroutine define_events

    !***************************************************************************
    !--- Allocate and assign the global alarm array
    !--- Associate alarms with events
    !***************************************************************************
    subroutine set_alarms(this_clock, start_time)

      !--- The clock to be associated with these alarms
      type(ESMF_Clock) :: this_clock

      !--- The start time for the clock
      type(ESMF_Time) :: start_time

      type(ESMF_Time)         :: alarm_time
      type(ESMF_TimeInterval) :: alarm_interval
      integer(ESMF_KIND_I)    :: interval_sec
      integer(ESMF_KIND_I)    :: rc
      integer(ESMF_KIND_I)    :: idx1
      character(64)  :: alarm_name = " "
      character(128) :: strng = " "

      !--- Coupler events must be defined prior to calling this routine
      if ( nevents <= 0 ) then
        write(6,*)"set_alarms: The events array is not defined."
        call xit("set_alarms",-1)
      endif

      if ( associated(alarm) ) nullify(alarm)
      allocate( alarm(nevents) )

      do idx1=1,nevents
        !--- Set the alarm name
        alarm_name = trim( event(idx1)%get_var )

        !--- Set the alarm interval (seconds)
        interval_sec = event(idx1)%freq

        if (verbose > 1) then
          write(6,'(a,i3,a,i8,2a)')' set_alarms: event ',idx1, &
              '  interval = ',interval_sec,'  action: ',trim(event(idx1)%action)
        endif

        call ESMF_TimeIntervalSet(alarm_interval, s=interval_sec, rc=rc)
        strng = "Problem setting alarm "//trim(alarm_name)//" time interval"
        call handle_ESMF_err(rc, trim(strng))

        !--- Set the alarm start time so that is does not ring until the first
        !--- alarm interval after the clock start time
        !xxx alarm_time = start_time + alarm_interval

        !--- Set the alarm start time
        alarm_time = start_time

        !--- Define the alarm
        alarm(idx1) = ESMF_AlarmCreate(clock=this_clock, ringTime=alarm_time, &
                       ringInterval=alarm_interval, name=trim(alarm_name), rc=rc)
        strng = "Problem creating alarm "//trim(alarm_name)
        call handle_ESMF_err(rc, trim(strng))

        !--- Associate the event alarm with the newly defined clock alarm
        event(idx1)%alarm_index = idx1
      enddo

    end subroutine set_alarms

    !***************************************************************************
    !--- Append an entry to regrid_list
    !***************************************************************************
    function append_regrid_list(new_regrid) result(new_id)
      type(regrid_t), pointer, intent(in) :: new_regrid
      integer :: new_id

      !--- Add a new regrid operation to the end of the list
      n_regrid_list = n_regrid_list + 1
      new_id = n_regrid_list
      if ( new_id > max_regrid_list ) then
        write(6,'(a)')"append_regrid_list: Maximum regrid list size exceeded."
        write(6,'(a)')"append_regrid_list: Increase max_regrid_list to ",new_id
        call xit("append_regrid_list",-1)
      endif

      !--- Assign the ID for this regrid operation
      new_regrid%regrid_id = new_id

      !--- Add the user supplied element to regrid_list(new_id)
      call copy_regrid_t( new_regrid, regrid_list(new_id) )

    end function append_regrid_list

    !***************************************************************************
    !--- Copy one regrid_t to another
    !***************************************************************************
    subroutine copy_regrid_t(regrid_in, regrid_out)
      type(regrid_t), intent(in)  :: regrid_in
      type(regrid_t), intent(out) :: regrid_out

      regrid_out % regrid_id       = regrid_in % regrid_id
      regrid_out % src_grid_id     = regrid_in % src_grid_id
      regrid_out % dst_grid_id     = regrid_in % dst_grid_id
      regrid_out % use_src_mask    = regrid_in % use_src_mask
      regrid_out % use_dst_mask    = regrid_in % use_dst_mask
      if ( associated(regrid_in%src_mask_values) ) then
        if ( associated(regrid_out%src_mask_values) ) deallocate(regrid_out%src_mask_values)
        allocate( regrid_out%src_mask_values(size(regrid_in%src_mask_values)) )
        regrid_out % src_mask_values = regrid_in % src_mask_values
      endif
      if ( associated(regrid_in%dst_mask_values) ) then
        if ( associated(regrid_out%dst_mask_values) ) deallocate(regrid_out%dst_mask_values)
        allocate( regrid_out%dst_mask_values(size(regrid_in%dst_mask_values)) )
        regrid_out % dst_mask_values = regrid_in % dst_mask_values
      endif
      regrid_out % regrid_method   = regrid_in % regrid_method
      regrid_out % norm_type       = regrid_in % norm_type
      regrid_out % line_type       = regrid_in % line_type
      regrid_out % src_grid_name   = regrid_in % src_grid_name
      regrid_out % dst_grid_name   = regrid_in % dst_grid_name
      regrid_out % src_grid_file_name = regrid_in % src_grid_file_name
      regrid_out % dst_grid_file_name = regrid_in % dst_grid_file_name
      regrid_out % wgt_grid_file_name = regrid_in % wgt_grid_file_name

      if ( associated(regrid_in%src_frac) ) then
        if ( associated(regrid_out%src_frac) ) deallocate(regrid_out%src_frac)
        allocate( regrid_out%src_frac(size(regrid_in%src_frac,dim=1), &
                                      size(regrid_in%src_frac,dim=2)) )
        regrid_out%src_frac = regrid_in%src_frac
      endif

      if ( associated(regrid_in%dst_frac) ) then
        if ( associated(regrid_out%dst_frac) ) deallocate(regrid_out%dst_frac)
        allocate( regrid_out%dst_frac(size(regrid_in%dst_frac,dim=1), &
                                      size(regrid_in%dst_frac,dim=2)) )
        regrid_out%dst_frac = regrid_in%dst_frac
      endif

      if ( associated(regrid_in%src_area) ) then
        if ( associated(regrid_out%src_area) ) deallocate(regrid_out%src_area)
        allocate( regrid_out%src_area(size(regrid_in%src_area,dim=1), &
                                      size(regrid_in%src_area,dim=2)) )
        regrid_out%src_area = regrid_in%src_area
      endif

      if ( associated(regrid_in%dst_area) ) then
        if ( associated(regrid_out%dst_area) ) deallocate(regrid_out%dst_area)
        allocate( regrid_out%dst_area(size(regrid_in%dst_area,dim=1), &
                                      size(regrid_in%dst_area,dim=2)) )
        regrid_out%dst_area = regrid_in%dst_area
      endif

      if ( associated(regrid_in%rh) ) then
        write(6,*)"copy_regrid_t: ????? Reassigning route handle. ????? "
        regrid_out%rh = regrid_in%rh
      else
        regrid_out%rh => null()
      endif
    end subroutine copy_regrid_t

    !***************************************************************************
    !--- Print elements of regrid_t
    !***************************************************************************
    subroutine print_regrid_t(rid)
      type(regrid_t), intent(in)  :: rid
      write(6,*)"            Regrid ID: ",rid % regrid_id
      write(6,*)"        Regrid Method: ",rid % regrid_method
      write(6,*)"            Norm Type: ",rid % norm_type
      write(6,*)"            Line Type: ",rid % line_type
      write(6,*)"       Source grid ID: ",rid % src_grid_id
      write(6,*)"         Dest grid ID: ",rid % dst_grid_id
      write(6,*)"      Use Source Mask: ",rid % use_src_mask
      if ( associated(rid%src_mask_values) ) then
        write(6,*)"   Source Mask Values: ",rid % src_mask_values
      endif
      write(6,*)"        Use Dest Mask: ",rid % use_dst_mask
      if ( associated(rid%dst_mask_values) ) then
        write(6,*)"     Dest Mask Values: ",rid % dst_mask_values
      endif
      write(6,*)"     Source grid name: ",rid % src_grid_name
      write(6,*)"       Dest grid name: ",rid % dst_grid_name
      write(6,*)"Source grid file name: ",rid % src_grid_file_name
      write(6,*)"  Dest grid file name: ",rid % dst_grid_file_name
      write(6,*)"    Weights file name: ",rid % wgt_grid_file_name
      if ( associated(rid%src_frac) ) then
        write(6,*)"Source fraction field is associated"
      else
        write(6,*)"Source fraction field is not associated"
      endif
      if ( associated(rid%dst_frac) ) then
        write(6,*)"Destination fraction field is associated"
      else
        write(6,*)"Destination fraction field is not associated"
      endif
      if ( associated(rid%src_area) ) then
        write(6,*)"Source area field is associated"
      else
        write(6,*)"Source area field is not associated"
      endif
      if ( associated(rid%dst_area) ) then
        write(6,*)"Destination area field is associated"
      else
        write(6,*)"Destination area field is not associated"
      endif
      if ( associated(rid%rh) ) then
        write(6,*)"Route Handle is associated"
      else
        write(6,*)"Route Handle is not associated"
      endif
      call flush(6)
    end subroutine print_regrid_t

    !***************************************************************************
    !--- Find the first regrid ID in regrid_list that matches the input
    !--- Optionally generate weights for a new regrid ID if an existing regrid
    !--- ID is not found
    !***************************************************************************
    subroutine find_regrid_id(id, src_grid_name, dst_grid_name, regrid_method, &
                 use_src_mask, use_dst_mask, src_mask_values, dst_mask_values, &
                 src_grid_id, dst_grid_id, norm_type, line_type, create)
      integer,      intent(out) :: id
      character(*), intent(in)  :: src_grid_name
      character(*), intent(in)  :: dst_grid_name
      character(*), intent(in)  :: regrid_method
      logical,      intent(in)  :: use_src_mask
      logical,      intent(in)  :: use_dst_mask
      integer,      intent(in)  :: src_mask_values(:)
      integer,      intent(in)  :: dst_mask_values(:)
      integer,      intent(in), optional :: src_grid_id
      integer,      intent(in), optional :: dst_grid_id
      character(*), intent(in), optional :: norm_type
      character(*), intent(in), optional :: line_type
      logical,      intent(in), optional :: create

      !--- Local
      integer :: idx1, idx2, verbose=3
      logical :: found, local_ignore_unmapped
      type(regrid_t), pointer :: new_regrid => null()
      real(kind=8), pointer :: cell_frac(:,:) => null()
      real(kind=8), pointer :: cpl_fland(:) => null()
      real(kind=8), pointer :: diff_fland(:) => null()
      real(kind=8), pointer :: alat(:) => null()
      real(kind=8), pointer :: alon(:) => null()
      real(kind=8) :: tol_fland=1e-6
      integer :: nxa, nya, ndiff
      integer(kind=4) :: ncfid, ncvid1, ncvid2, ncvid3
      integer(kind=4) :: ncxid, ncx_dim_id, ncyid, ncy_dim_id

      !--- Return zero if no matches are found
      id = 0

      !--- Ensure input character variables are not blank
      if ( len_trim(src_grid_name) < 1 ) then
        write(6,*)"find_regrid_id: src_grid_name is blank"
        call flush(6)
        call xit("find_regrid_id",-1)
      endif
      if ( len_trim(dst_grid_name) < 1 ) then
        write(6,*)"find_regrid_id: dst_grid_name is blank"
        call flush(6)
        call xit("find_regrid_id",-2)
      endif
      if ( len_trim(regrid_method) < 1 ) then
        write(6,*)"find_regrid_id: regrid_method is blank"
        call flush(6)
        call xit("find_regrid_id",-3)
      endif
      if ( present(norm_type) ) then
        if ( len_trim(norm_type) < 1 ) then
          write(6,*)"find_regrid_id: norm_type is blank"
          call flush(6)
          call xit("find_regrid_id",-4)
        endif
      endif
      if ( present(line_type) ) then
        if ( len_trim(line_type) < 1 ) then
          write(6,*)"find_regrid_id: line_type is blank"
          call flush(6)
          call xit("find_regrid_id",-5)
        endif
      endif

      if ( present(src_grid_id) ) then
        !--- Verify that the input src grid ID is within the valid range
        if ( src_grid_id < 1 .or. src_grid_id > n_grid_list ) then
          write(6,*)"find_regrid_id: Invalid src_grid_id = ",src_grid_id
          call flush(6)
          call xit("find_regrid_id",-1)
        endif
      endif

      if ( present(dst_grid_id) ) then
        !--- Verify that the input dst grid ID is within the valid range
        if ( dst_grid_id < 1 .or. dst_grid_id > n_grid_list ) then
          write(6,*)"find_regrid_id: Invalid dst_grid_id = ",dst_grid_id
          call flush(6)
          call xit("find_regrid_id",-2)
        endif
      endif

      if ( n_regrid_list > 0 ) then
        !--- At least one regrid operation is already defined
        !--- Search for an existing regrid ID that matches the input request
        do idx1=1,n_regrid_list
          found = .true.

          !--- Compare input src grid with that in the regrid list
          if ( present(src_grid_id) ) then
            !--- If src grid ID is present then compare it and ignore src_grid_name
            if ( regrid_list(idx1)%src_grid_id .ne. src_grid_id ) found = .false.
          else
            !--- Otherwise use the grid name to determine grid ID below
            if ( trim(adjustl(src_grid_name)) .ne. &
                 trim(adjustl(regrid_list(idx1)%src_grid_name)) ) found = .false.
          endif

          !--- Compare input dst grid with that in the regrid list
          if ( present(dst_grid_id) ) then
            !--- If dst grid ID is present then compare it and ignore dst_grid_name
            if ( regrid_list(idx1)%dst_grid_id .ne. dst_grid_id ) found = .false.
          else
            !--- Otherwise use the grid name to determine grid ID below
            if ( trim(adjustl(dst_grid_name)) .ne. &
                 trim(adjustl(regrid_list(idx1)%dst_grid_name)) ) found = .false.
          endif

          !--- Compare input regrid method with that in the regrid list
          if ( trim(adjustl(regrid_method)) .ne. &
               trim(adjustl(regrid_list(idx1)%regrid_method)) )  found=.false.
          if ( present(norm_type) ) then
            if ( trim(adjustl(norm_type)) .ne. &
                 trim(adjustl(regrid_list(idx1)%norm_type)) ) found=.false.
          endif

          if ( present(line_type) ) then
            if ( trim(adjustl(line_type)) .ne. &
                 trim(adjustl(regrid_list(idx1)%line_type)) ) found=.false.
          endif

          !--- Compare input src mask options with that in the regrid list
          if ( use_src_mask .neqv. regrid_list(idx1)%use_src_mask ) found=.false.
          if ( use_src_mask .and. associated(regrid_list(idx1)%src_mask_values) ) then
            if ( size(src_mask_values) == size(regrid_list(idx1)%src_mask_values) ) then
              if ( any(src_mask_values /= regrid_list(idx1)%src_mask_values) ) found=.false.
            else
              found = .false.
            endif
          endif

          !--- Compare input dst mask options with that in the regrid list
          if ( use_dst_mask .neqv. regrid_list(idx1)%use_dst_mask ) found=.false.
          if ( use_dst_mask .and. associated(regrid_list(idx1)%dst_mask_values) ) then
            if ( size(dst_mask_values) == size(regrid_list(idx1)%dst_mask_values) ) then
              if ( any(dst_mask_values /= regrid_list(idx1)%dst_mask_values) ) found=.false.
            else
              found = .false.
            endif
          endif

          if ( found) then
            !--- Accept the first match found
            id = regrid_list(idx1)%regrid_id
            exit
          endif
        enddo
      endif

      if ( present(create) ) then
        !--- If the regrid operation was not found then define it when create = T

        if ( id == 0 .and. create) then
          !--- Define a new regrid operation
          allocate( new_regrid )

          if ( present(src_grid_id) ) then
            new_regrid%src_grid_name = trim(cpl_grid(src_grid_id)%file_name)
          else
            new_regrid%src_grid_name = trim(src_grid_name)
          endif

          if ( present(dst_grid_id) ) then
            new_regrid%dst_grid_name = trim(cpl_grid(dst_grid_id)%file_name)
          else
            new_regrid%dst_grid_name = trim(dst_grid_name)
          endif

          new_regrid%regrid_method = trim(regrid_method)

          if ( present(norm_type) ) then
            new_regrid%norm_type = trim(norm_type)
          else
            new_regrid%norm_type = "fracarea"
          endif
          if ( present(line_type) ) then
            new_regrid%line_type = trim(line_type)
          else
            new_regrid%line_type = "spherical"
          endif

          !--- Identify the src grid Id
          if ( present(src_grid_id) ) then
            new_regrid%src_grid_id = src_grid_id
          else
            new_regrid%src_grid_id = find_grid_id( src_grid_name )
          endif
          if ( new_regrid%src_grid_id == 0 ) then
            write(6,*)"find_regrid_id: Undefined src grid. name = ",trim(src_grid_name)
            call flush(6)
            call xit("find_regrid_id",-3)
          endif

          !--- Identify the dst grid Id
          if ( present(dst_grid_id) ) then
            new_regrid%dst_grid_id = dst_grid_id
          else
            new_regrid%dst_grid_id = find_grid_id( dst_grid_name )
          endif
          if ( new_regrid%dst_grid_id == 0 ) then
            write(6,*)"find_regrid_id: Undefined dst grid. name = ",trim(dst_grid_name)
            call flush(6)
            call xit("find_regrid_id",-4)
          endif

          !--- Identify the src grid file name
          new_regrid%src_grid_file_name = trim( grid_desc_file_name(new_regrid%src_grid_id) )

          !--- Identify the dst grid file name
          new_regrid%dst_grid_file_name = trim( grid_desc_file_name(new_regrid%dst_grid_id) )

          new_regrid%use_src_mask = use_src_mask
          new_regrid%use_dst_mask = use_dst_mask

          allocate( new_regrid%src_mask_values(size(src_mask_values)) )
          new_regrid%src_mask_values = src_mask_values
          allocate( new_regrid%dst_mask_values(size(dst_mask_values)) )
          new_regrid%dst_mask_values = dst_mask_values

          !--- Add this to the regrid list and return the new id
          id = append_regrid_list(new_regrid)

          nullify( new_regrid )

          if ( id == 0 ) then
            !--- If id is still zero then we have a problem
            write(6,*)'find_regrid_id: Unable to define a new regid operation.'
            call flush(6)
            call xit("find_regrid_id",-5)
          endif

          if ( verbose > 2 ) then
            write(6,*)"find_regrid_id: New regrid ID=",id
            call print_regrid_t( regrid_list(id) )
            call flush(6)
          endif

          local_ignore_unmapped = .true.

          !--- Define the new regrid operation
          call define_regrid(id, ignore_unmapped=local_ignore_unmapped )

          if ( verbose > 1 ) then
            write(6,'(7a,i3)') &
                "find_regrid_id: New regrid operation is defined for a ", &
                trim(regrid_method),' remap from grid ', &
                trim(src_grid_name)," to grid ",trim(dst_grid_name),"  ID=",id
            call flush(6)
          endif

          !--- Check certain features of this remap
          call check_remap_global_sum(id)

          if ( regrid_list(id)%src_grid_id == id_atm_grid .and. &
               regrid_list(id)%dst_grid_id == id_ocn_grid ) then
            !--- This is a remap operation from the AGCM to the OCEAN

            !--- Extract the src_frac from this regrid operation to compare with fland
            if ( associated(atm_fland) ) then
              if ( associated(regrid_list(id)%src_frac) ) then
                !--- Assign cell_frac from the src area fraction
                nxa = size(regrid_list(id)%src_frac,dim=1)
                nya = size(regrid_list(id)%src_frac,dim=2)

                !--- Add 1 overlap longitude to be consistent with atm_fland
                if ( associated(cell_frac) ) deallocate(cell_frac)
                allocate( cell_frac(nxa+1,nya) )
                cell_frac(1:nxa,:) = regrid_list(id)%src_frac(:,:)
                cell_frac(nxa+1,:) = cell_frac(1,:)

                !--- Reshape cell_frac as a 1D array, like atm_fland
                if ( associated(cpl_fland) ) deallocate(cpl_fland)
                allocate( cpl_fland(size(cell_frac)) )
                cpl_fland = reshape( cell_frac, (/size(cell_frac)/) )

                !--- src_frac may contain values very close to 1
                !--- These cells will be entirely ocean so we set them
                !--- to be exactly 1 for use in fland
                where( cpl_fland > 1.0_8-tol_fland ) cpl_fland = 1.0_8

                !--- src_frac may contain values very close to 0
                !--- These cells will be entirely land so we set them
                !--- to be exactly 0 for use in fland
                where( cpl_fland < tol_fland ) cpl_fland = 0.0_8

                !--- Invert land/ocean values so that 1=land and 0=ocean
                cpl_fland = 1.0_8 - cpl_fland
                where( cpl_fland > 1.0_8-tol_fland ) cpl_fland = 1.0_8
                where( cpl_fland < tol_fland )       cpl_fland = 0.0_8

                if ( size(cpl_fland) == size(atm_fland) ) then
                  !--- Write fland to a local file
                  !--- Put the difference into cpl_fland
                  if ( associated(diff_fland) ) deallocate(diff_fland)
                  allocate( diff_fland(size(atm_fland)) )
                  diff_fland = atm_fland - cpl_fland
                  where (diff_fland < 0.0_8) diff_fland = -diff_fland  ! abs(diff_fland)

                  !--- Count the number of grid cells the have different fland values between
                  !--- those passed from the AGCM and those calculated in the coupler
                  ndiff = count(diff_fland > tol_fland)
                  write(6,*)"find_regrid_id: Max of abs(atm_fland-cpl_fland) is ",maxval(diff_fland)

                  if ( ndiff > 0 ) then
                    !--- If any AGCM cells have a different fland value between that which
                    !--- is passed from the AGCM and that which is calculated in the coupler
                    !--- then flag as a fatal error
                    write(6,*)"find_regrid_id: **EE** ",ndiff," FLAND values are inconsistent with NEMO tmask"
                    write(6,*)"    FLAND maximum difference  ",maxval(diff_fland)
                    !--- Recycle cell_frac as tmp space
                    cell_frac = 0.0_8
                    cell_frac = reshape( diff_fland, (/nxa+1,nya/) )
                    do idx1=1,size(cell_frac,dim=1)
                      do idx2=1,size(cell_frac,dim=2)
                        if ( abs(cell_frac(idx1,idx2)) > tol_fland ) then
                          if ( idx1 == nxa+1 ) then
                            write(6,'(2i6,a,2f12.6,a,g16.8)')idx1,idx2, &
                                "  lon,lat=",atm_lon(1),atm_lat(idx2), &
                                "  abs(diff)=",cell_frac(idx1,idx2)
                          else
                            write(6,'(2i6,a,2f12.6,a,g16.8)')idx1,idx2, &
                                "  lon,lat=",atm_lon(idx1),atm_lat(idx2), &
                                "  abs(diff)=",cell_frac(idx1,idx2)
                          endif
                        endif
                      enddo
                    enddo
                    call flush(6)

                    !--- Dump the atm_fland, cpl_fland and their difference to a netcdf file
                    call ncdf_open(ncfid, "cpl_fland_info.nc", "create")
                    call get_atm_lon_1d(alon, overlap=.true.)
                    call get_atm_lat_1d(alat)
                    call ncdf_add_coord(ncfid, ncxid, ncx_dim_id, "lon", vals=alon)
                    call ncdf_add_coord(ncfid, ncyid, ncy_dim_id, "lat", vals=alat)
                    cell_frac = 0.0_8
                    cell_frac = reshape(atm_fland, (/nxa+1,nya/))
                    call ncdf_add_var(ncfid, ncvid1, "fland_atm", (/ncxid,ncyid/), vals=cell_frac)
                    call ncdf_add_att(ncfid, ncvid1, "long_name", "fland received from the AGCM")
                    cell_frac = reshape(cpl_fland, (/nxa+1,nya/))
                    call ncdf_add_var(ncfid, ncvid2, "fland_cpl", (/ncxid,ncyid/), vals=cell_frac)
                    call ncdf_add_att(ncfid, ncvid2, "long_name", "fland calculated in the coupler")
                    cell_frac = reshape(diff_fland, (/nxa+1,nya/))
                    call ncdf_add_var(ncfid, ncvid3, "fland_diff", (/ncxid,ncyid/), vals=cell_frac)
                    call ncdf_add_att(ncfid, ncvid3, "long_name", "difference fland_atm - fland_cpl")
                    call ncdf_close(ncfid)
                    write(6,*)"find_regrid_id: Created file cpl_fland_info.nc containing FLAND values"

                    !--- Dump the fland value that is calculated in the coupler to a
                    !--- CCCma binary file (with values inverted so that 1=ocean, 0=land)
                    cpl_fland = 1.0_8 - cpl_fland
                    where( cpl_fland > 1.0_8 ) cpl_fland = 1.0_8
                    where( cpl_fland < 1e-14 ) cpl_fland = 0.0_8
                    call dump_farray("SRCF", cpl_fland, nx=nxa+1, ny=nya, fname="cpl_fland")
                    write(6,*)"find_regrid_id: Created CCCma binary file cpl_fland"
                    write(6,*)"  This file contains the value of fland calculated in the coupler"
                    call flush(6)
                    if (landfrac_bug) call xit("find_regrid_id",-6)
                  endif
                else
                  write(6,*)"find_regrid_id: Size mismatch between src_frac and atm_fland"
                  write(6,*)"    expected ",size(atm_fland)
                  write(6,*)"       found ",size(cpl_fland)
                  write(6,*)"find_regrid_id: ** WW ** Cannot validate AGCM fractional land mask"
                  call flush(6)
                endif
              else
                write(6,*)"find_regrid_id: src_frac is not allocated for ATM to OCN regrid operation"
                write(6,*)"find_regrid_id: ** WW ** Cannot validate AGCM fractional land mask"
                call flush(6)
              endif
            else
              write(6,*)"find_regrid_id: atm_fland is not allocated"
              write(6,*)"find_regrid_id: ** WW ** Cannot validate AGCM fractional land mask"
              call flush(6)
            endif

            if ( associated(cell_frac) )  deallocate(cell_frac)
            if ( associated(cpl_fland) )  deallocate(cpl_fland)
            if ( associated(diff_fland) ) deallocate(diff_fland)
          endif

        endif
      endif

      if ( id == 0 .and. verbose > 0 ) then
        write(6,'(6a)')" **WW** find_regrid_id: No regridding operation is defined for a ", &
            trim(regrid_method),' remap from grid ', &
            trim(src_grid_name)," to grid ",trim(dst_grid_name)
        call flush(6)
      endif

    end subroutine find_regrid_id

    !***************************************************************************
    !--- Evaluate source and destination global areas for this regrid operation
    !***************************************************************************
    subroutine check_remap_global_sum(regridID)
      integer, intent(in) :: regridID

      !--- Local
      integer :: nxa, nya, src_grid_id, dst_grid_id
      real(kind=8), pointer :: src_mask(:,:) => null(),  dst_mask(:,:) => null()
      real(kind=8), pointer :: src_area(:,:) => null(),  dst_area(:,:) => null()
      real(kind=8), pointer :: src_frac(:,:) => null(),  dst_frac(:,:) => null()
      real(kind=8), pointer :: src_area_esmf(:,:) => null(), dst_area_esmf(:,:) => null()
      real(kind=8), pointer :: awrk2d(:,:) => null(), bwrk2d(:,:) => null()
      real(kind=8), parameter :: pi4 = 12.56637061435917_8
      character(32) :: src_grid_name, dst_grid_name

      nullify(src_mask, dst_mask, src_area, dst_area, src_frac, dst_frac)
      nullify(src_area_esmf, dst_area_esmf, awrk2d, bwrk2d)

      if ( regridID < 1 .or. regridID > n_regrid_list ) then
        write(6,'(a,i3)')"check_remap_global_sum: Regrid ID is out of range. ID=",regridID
        call xit("check_remap_global_sum",-1)
      endif

      if ( regrid_list(regridID)%src_grid_id == id_atm_grid .and. &
           regrid_list(regridID)%dst_grid_id == id_ocn_grid ) then
        write(6,*)"check_remap_global_sum: Remap from AGCM grid to NEMO grid"
        src_grid_name = "AGCM"
        dst_grid_name = "NEMO"
      else if ( regrid_list(regridID)%src_grid_id == id_ocn_grid .and. &
                regrid_list(regridID)%dst_grid_id == id_atm_grid ) then
        write(6,*)"check_remap_global_sum: Remap from NEMO grid to AGCM grid"
        src_grid_name = "NEMO"
        dst_grid_name = "AGCM"
      else
        write(6,*)"check_remap_global_sum: Remap from ", &
            trim(regrid_list(regridID)%src_grid_name)," grid to ", &
            trim(regrid_list(regridID)%dst_grid_name)," grid"
        src_grid_name = trim(regrid_list(regridID)%src_grid_name)
        dst_grid_name = trim(regrid_list(regridID)%dst_grid_name)
      endif

      if ( associated(regrid_list(regridID)%src_frac) .and. &
           associated(regrid_list(regridID)%src_area) ) then
        !--- Assign weights for the global integral from the src grid
        nxa = size(regrid_list(regridID)%src_frac,dim=1)
        nya = size(regrid_list(regridID)%src_frac,dim=2)
        allocate( src_frac(nxa,nya), src_area(nxa,nya) )
        src_frac = regrid_list(regridID)%src_frac
        src_area = regrid_list(regridID)%src_area
      else
        write(6,*)"check_remap_global_sum: src_frac or src_area is missing"
        call flush(6)
        return
      endif

      if ( associated(regrid_list(regridID)%dst_frac) .and. &
           associated(regrid_list(regridID)%dst_area) ) then
        !--- Assign weights for the global integral from the dst grid
        nxa = size(regrid_list(regridID)%dst_frac,dim=1)
        nya = size(regrid_list(regridID)%dst_frac,dim=2)
        allocate( dst_frac(nxa,nya), dst_area(nxa,nya) )
        dst_frac = regrid_list(regridID)%dst_frac
        dst_area = regrid_list(regridID)%dst_area
      else
        write(6,*)"check_remap_global_sum: dst_frac or dst_area is missing"
        call flush(6)
        return
      endif

      !--- Identify the source and destination grids by ID
      src_grid_id = regrid_list(regridID)%src_grid_id
      dst_grid_id = regrid_list(regridID)%dst_grid_id

      !--- Retrieve the grid masks associated with src and dst grids
      call get_grid_item(src_mask, src_grid_id, "center", "mask")
      call get_grid_item(dst_mask, dst_grid_id, "center", "mask")

      !--- Calculate the grid cell area via a call to ESMF_FieldRegridGetArea
      call compute_grid_area(src_grid_id, src_area_esmf)
      call compute_grid_area(dst_grid_id, dst_area_esmf)

      !--- If the data field is identically 1 then the source and destination
      !--- global integrals are as follows
      write(6,*)"check_remap_global_sum: Global integral of unit function on ",trim(src_grid_name)," grid (frac*area)"
      write(6,*)"check_remap_global_sum:     with areas that are used in remap                = ",sum(src_frac*src_area)
      write(6,*)"check_remap_global_sum:     with areas calculated by ESMF_FieldRegridGetArea = ",sum(src_frac*src_area_esmf)
      write(6,*)"check_remap_global_sum: Global integral of unit function on ",trim(dst_grid_name)," grid (frac*area)"
      write(6,*)"check_remap_global_sum:     with areas that are used in remap                = ",sum(dst_frac*dst_area)
      write(6,*)"check_remap_global_sum:     with areas calculated by ESMF_FieldRegridGetArea = ",sum(dst_frac*dst_area_esmf)

      write(6,*)"check_remap_global_sum: global integral of src mask*area  = ",sum(src_mask*src_area)
      write(6,*)"check_remap_global_sum: global integral of dst mask*area  = ",sum(dst_mask*dst_area)
      write(6,*)"check_remap_global_sum:              area of unit sphere  = ",pi4

      if ( associated(awrk2d) ) deallocate(awrk2d)
      allocate( awrk2d(size(src_mask,dim=1),size(src_mask,dim=2)) )
      awrk2d = src_mask-src_frac
      write(6,*)"check_remap_global_sum: min(src_mask-src_frac)  = ",minval(awrk2d)
      write(6,*)"check_remap_global_sum: max(src_mask-src_frac)  = ",maxval(awrk2d)
      write(6,*)"check_remap_global_sum: avg(src_mask-src_frac)  = ",sum(awrk2d)/real(size(awrk2d),kind=8)

      if ( associated(awrk2d) ) deallocate(awrk2d)
      allocate( awrk2d(size(dst_mask,dim=1),size(dst_mask,dim=2)) )
      awrk2d = dst_mask-dst_frac
      write(6,*)"check_remap_global_sum: min(dst_mask-dst_frac)  = ",minval(awrk2d)
      write(6,*)"check_remap_global_sum: max(dst_mask-dst_frac)  = ",maxval(awrk2d)
      write(6,*)"check_remap_global_sum: avg(dst_mask-dst_frac)  = ",sum(awrk2d)/real(size(awrk2d),kind=8)
      call flush(6)

      if ( associated(src_mask) )  deallocate(src_mask)
      if ( associated(dst_mask) )  deallocate(dst_mask)
      if ( associated(src_area) )  deallocate(src_area)
      if ( associated(dst_area) )  deallocate(dst_area)
      if ( associated(src_frac) )  deallocate(src_frac)
      if ( associated(dst_frac) )  deallocate(dst_frac)
      if ( associated(awrk2d) )    deallocate(awrk2d)
      if ( associated(bwrk2d) )    deallocate(bwrk2d)

    end subroutine check_remap_global_sum

    !***************************************************************************
    !--- Create a route handle for a particular regrid operation
    !--- This route handle will be associated with the id_regrid index of the
    !--- global array regrid_rh_list(:)
    !***************************************************************************
    subroutine define_regrid(id_regrid, use_src_mask, use_dst_mask, ignore_unmapped)
#undef save_regrid_weights
#ifdef save_regrid_weights
      use ESMF_RegridWeightGenMod
      use ESMF_RegridWeightGenCheckMod
#endif
      integer, intent(in) :: id_regrid
      logical, intent(in), optional :: use_src_mask
      logical, intent(in), optional :: use_dst_mask
      logical, intent(in), optional :: ignore_unmapped

      !--- Local
      integer :: id_src_grid, id_dst_grid
      integer :: grid_file_type
      character(128) :: map_method, src_grid_name, dst_grid_name, regrid_line_type
      type(ESMF_grid)  :: src_grid, dst_grid
      type(ESMF_Field) :: src_field, dst_field
      type(ESMF_Field) :: src_FracField, dst_FracField

      type(ESMF_RegridMethod_Flag)   ::   regrid_method
      type(ESMF_NormType_Flag)       ::       norm_type
      type(ESMF_PoleMethod_Flag)     ::     pole_method
      integer(ESMF_KIND_I)           ::       pole_npts = 0
      type(ESMF_UnmappedAction_Flag) :: unmapped_action
      type(ESMF_LineType_Flag)       ::       line_type

      real   (ESMF_KIND_R8), pointer, dimension(:)   :: factor_list => null()
      integer(ESMF_KIND_I4), pointer, dimension(:,:) :: factor_index_list => null()

      integer(ESMF_KIND_I4), pointer, dimension(:) :: unmapped_dst_list => null()

      !--- Pointer used to assign local arrays
      real(ESMF_KIND_R8), pointer, dimension(:,:) :: farray_ptr => null()
      real(ESMF_KIND_R8), pointer, dimension(:,:) :: src_area => null()
      real(ESMF_KIND_R8), pointer, dimension(:,:) :: dst_area => null()

      !--- Pointers used to assign local coordinate arrays
      real(ESMF_KIND_R8), pointer, dimension(:,:) :: lon_c => null()
      real(ESMF_KIND_R8), pointer, dimension(:,:) :: lat_c => null()
      real(ESMF_KIND_R8), pointer, dimension(:,:) :: lon_v => null()
      real(ESMF_KIND_R8), pointer, dimension(:,:) :: lat_v => null()
      integer :: nx_c, ny_c, nx_v, ny_v, ix, iy, ixp1, iyp1

      real(ESMF_KIND_R8), pointer :: coord2d(:,:) => null()
      real(ESMF_KIND_R8), pointer :: coord1d(:) => null()
      type(ESMF_Array) :: coordAry
      integer(ESMF_KIND_I) :: excCount(2)

      !--- 1D version of grid area
      real(ESMF_KIND_R8), pointer, dimension(:) :: srcArea_out => null()
      real(ESMF_KIND_R8), pointer, dimension(:) :: dstArea_out => null()

      !--- 1D version of grid fraction
      real(ESMF_KIND_R8), pointer, dimension(:) :: srcFrac_out => null()
      real(ESMF_KIND_R8), pointer, dimension(:) :: dstFrac_out => null()

      !--- Pointer used to access data in field objects
      real(ESMF_KIND_R8), pointer, dimension(:,:) :: src_data_ptr => null()

      integer(ESMF_KIND_I) :: coord_dim, rank
      integer(ESMF_KIND_I), pointer, dimension(:) :: src_mask_values => null()
      integer(ESMF_KIND_I), pointer, dimension(:) :: dst_mask_values => null()

      integer(ESMF_KIND_I) :: src_term_processing = 0
      integer(ESMF_KIND_I) :: pipeline_depth = 1

      integer(ESMF_KIND_I) :: rc, lnum
      character(64) :: method = "define_regrid"
      logical :: local_use_src_mask, local_use_dst_mask
      logical :: local_ignore_unmapped
      logical :: exists
      logical :: save_weights_to_file
      integer :: idx1, idx2, n_unmapped
      integer :: nx, ny, nx_src, ny_src, nx_dst, ny_dst
      character(256) :: wgt_file, wgt_title, this_norm_type
      character(256) :: src_grid_file_name, dst_grid_file_name
      character(256) :: src_grid_title, dst_grid_title
      real :: time_in, time_out
      type(cpl_vinfo_t), pointer :: tmp_vinfo => null()
      character(128) :: strng = " "

      !--- Determine if a mask will be used on the source grid
      !--- Masks are read from the grid description file
      if ( present(use_src_mask) ) then
        local_use_src_mask = use_src_mask
      else
        local_use_src_mask = regrid_list(id_regrid)%use_src_mask
      endif

      !--- Determine if a mask will be used on the destination grid
      !--- Masks are read from the grid description file
      if ( present(use_dst_mask) ) then
        local_use_dst_mask = use_dst_mask
      else
        local_use_dst_mask = regrid_list(id_regrid)%use_dst_mask
      endif

      !--- Should we ignore umapped values or not
      if ( present(ignore_unmapped) ) then
        local_ignore_unmapped = ignore_unmapped
      else
        local_ignore_unmapped = .false.
      endif

      if ( id_regrid < 1 .or. id_regrid > n_regrid_list ) then
        write(6,'(a,i3)')"define_regrid: Regrid ID is out of range. ID=",id_regrid
        call xit("define_regrid",-1)
      endif

      !--- Assign local grid IDs and regrid method from regrid_list
      id_src_grid = regrid_list(id_regrid)%src_grid_id
      id_dst_grid = regrid_list(id_regrid)%dst_grid_id
      map_method = trim( regrid_list(id_regrid)%regrid_method )
      this_norm_type = trim( regrid_list(id_regrid)%norm_type )
      regrid_line_type = trim( regrid_list(id_regrid)%line_type )

      !--- Assign local src and dst grid objects from the global list
      src_grid = grid_list(id_src_grid)
      dst_grid = grid_list(id_dst_grid)

      !--- Assign local src and dst grid names
      src_grid_name = trim( regrid_list(id_regrid)%src_grid_name )
      dst_grid_name = trim( regrid_list(id_regrid)%dst_grid_name )

      !--- Assign local src and dst grid file names
      src_grid_file_name = trim( regrid_list(id_regrid)%src_grid_file_name )
      dst_grid_file_name = trim( regrid_list(id_regrid)%dst_grid_file_name )

      write(6,'(4a)')'define_regrid: SRC grid: ',trim(src_grid_name), &
          '  file name: ',trim(src_grid_file_name)
      write(6,'(5a)')'define_regrid: DST grid: ',trim(dst_grid_name), &
          '  file name: ',trim(dst_grid_file_name)
      write(6,'(2a)')'define_regrid: Regrid method: ',trim(map_method)
      write(6,'(2a)')'define_regrid: Norm type: ',trim(this_norm_type)

      !--- Create ESMF field objects from an ESMF grid and typekind

      !--- src_field will hold the source field values
      src_field = ESMF_FieldCreate( grid_list(id_src_grid), ESMF_TYPEKIND_R8, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Problem creating src field", line=lnum, method=method)

      !--- The source field MUST contain valid data (ie internal data must be initialized)
      !--- If it does not contain valid data the call to ESMF_FieldRegridStore below
      !--- will result in a segfault
      if ( associated(src_data_ptr) ) nullify( src_data_ptr )
      call ESMF_FieldGet( src_field, farrayPtr=src_data_ptr, rc=rc )
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Problem getting src_field pointer", line=lnum, method=method)
      if ( associated(src_data_ptr) ) then
        !--- Initialize internal data to zero to avoid a segfault
        src_data_ptr = 0.0
        nullify(src_data_ptr)
      endif

      !--- src_fracfield will hold the fraction of each source cell
      !--- participating in the regridding
      src_fracfield = ESMF_FieldCreate( src_grid, ESMF_TYPEKIND_R8, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Problem creating src fracfield", line=lnum, method=method)

      !--- dst_field will hold the destination field values
      dst_field = ESMF_FieldCreate( grid_list(id_dst_grid), ESMF_TYPEKIND_R8, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Problem creating dst field", line=lnum, method=method)

      !--- dst_fracfield will hold the fraction of each destination cell
      !--- participating in the regridding
      dst_fracfield = ESMF_FieldCreate( dst_grid, ESMF_TYPEKIND_R8, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Problem creating dst fracfield", line=lnum, method=method)

      if ( verbose > 10 ) then
        call ESMF_FieldPrint(src_field, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Problem printing src field", line=lnum, method=method)

        call ESMF_FieldPrint(dst_field, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Problem printing dst field", line=lnum, method=method)
      endif

      !--- Valid regrid methods are
      !--- ESMF_REGRIDMETHOD_BILINEAR     ...Bilinear interpolation
      !--- ESMF_REGRIDMETHOD_PATCH        ...Higher-order patch recovery interpolation
      !--- ESMF_REGRIDMETHOD_NEAREST_STOD ...Each destination point is mapped
      !---                                   to the closest source point
      !--- ESMF_REGRIDMETHOD_NEAREST_DTOS ...Each source point is mapped to the
      !---                                   closest destination point
      !--- ESMF_REGRIDMETHOD_CONSERVE     ...First order conservative interpolation

      !--- Assign values to regrid_method, pole_method and pole_npts
      select case( trim(map_method) )
      case ('conservative', 'conserve', 'cv')
        regrid_method = ESMF_REGRIDMETHOD_CONSERVE
        pole_method = ESMF_POLEMETHOD_NONE
        pole_npts   = 0
      case ('bilinear', 'bl')
        regrid_method = ESMF_REGRIDMETHOD_BILINEAR
        pole_method = ESMF_POLEMETHOD_ALLAVG
        pole_npts   = -1
      case ('patch')
        regrid_method = ESMF_REGRIDMETHOD_PATCH
        pole_method = ESMF_POLEMETHOD_ALLAVG
        pole_npts   = -1
      case ('nearest_stod')
        regrid_method = ESMF_REGRIDMETHOD_NEAREST_STOD
        pole_method = ESMF_POLEMETHOD_ALLAVG
        pole_npts   = -1
      case ('nearest_dtos')
        regrid_method = ESMF_REGRIDMETHOD_NEAREST_DTOS
        pole_method = ESMF_POLEMETHOD_ALLAVG
        pole_npts   = -1
      case default
        write(6,*)'define_regrid: Unknown mapping method ',trim(map_method)
        call xit("define_regrid",-1)
      end select

      !--- line_type allows the user to select the path of the line which connects two
      !--- points on the surface of a sphere. This in turn controls the path along which
      !--- distances are calculated and the shape of the edges that make up a cell
      !--- ESMF_LINETYPE_CART         ...Cartesian line
      !--- ESMF_LINETYPE_GREAT_CIRCLE ...Great circle line

      !--- Define line_type
      select case( trim(regrid_line_type) )
      case ('spherical')
        line_type = ESMF_LINETYPE_GREAT_CIRCLE
      case ('cartesian')
        line_type = ESMF_LINETYPE_CART
        !--- Always use ESMF_POLEMETHOD_NONE on a cartesian grid
        pole_method = ESMF_POLEMETHOD_NONE
      case default
        write(6,*)'define_regrid: Unknown line tpe '
        call xit("define_regrid",-1)
      end select

      !--- For conservative interpolation the normalization options are
      !--- ESMF_NORMTYPE_DSTAREA  ...Destination area normalization
      !--- ESMF_NORMTYPE_FRACAREA ...Fraction area normalization

      !--- Assign a value to norm_type
      select case( trim(this_norm_type) )
      case ('fracarea')
        norm_type = ESMF_NORMTYPE_FRACAREA
      case ('destarea')
        norm_type = ESMF_NORMTYPE_DSTAREA
      case default
        write(6,*)'define_regrid: Unknown normalization option ',trim(this_norm_type)
        call xit("define_regrid",-1)
      end select

      !--- Possible actions taken when a destination point is unpmapped are
      !--- ESMF_UNMAPPEDACTION_ERROR  ...An error is issued when there exist destination
      !---                               points in a regridding operation that are not
      !---                               mapped by corresponding source points
      !--- ESMF_UNMAPPEDACTION_IGNORE ...Destination points which do not have corresponding
      !---                               source points are ignored and zeros are used for
      !---                               the entries of the sparse matrix that is generated

      !--- Decide what to do when a destination point cannot be mapped to any
      !--- source grid cell
      if (local_ignore_unmapped) then
        !--- Ignore these points
        unmapped_action = ESMF_UNMAPPEDACTION_IGNORE
      else
        !--- Abort if any of these points are found
        unmapped_action = ESMF_UNMAPPEDACTION_ERROR
      endif

      call cpu_time( time_in )

      !--- A location is masked if and only if the value for that location in the mask
      !--- that is read from the grid description file matches one of the values listed
      !--- in srcMaskValues or dstMaskValues
      !--- The absence of either srcMaskValues or dstMaskValues indicates that
      !--- no masking should be done for that Field ie
      !---   no srcMaskValue parameter indicates that source masking will not occur
      !---   no dstMaskValue parameter indicates that destination masking will not occur
      !--- It is assumed in the following that a value of 0 in imask, read from the grid
      !--- description file, will indicate that the point is masked

      !--- Assign mask values to use for src and dst fields (normally zero)

      nullify(src_mask_values)
      if ( associated(regrid_list(id_regrid)%src_mask_values) ) then
        if ( associated(src_mask_values) ) deallocate(src_mask_values)
        allocate( src_mask_values(size(regrid_list(id_regrid)%src_mask_values)) )
        src_mask_values = regrid_list(id_regrid)%src_mask_values
      endif
      nullify(dst_mask_values)
      if ( associated(regrid_list(id_regrid)%dst_mask_values) ) then
        if ( associated(dst_mask_values) ) deallocate(dst_mask_values)
        allocate( dst_mask_values(size(regrid_list(id_regrid)%dst_mask_values)) )
        dst_mask_values = regrid_list(id_regrid)%dst_mask_values
      endif

      !--- Create the weights and store route handle in regrid_rh_list
      if ( local_use_src_mask .and. local_use_dst_mask ) then
        !--- Mask both source and destination fields
        if ( .not. associated(src_mask_values) ) then
          write(6,*)'define_regrid: Masking src field but src_mask_values is not allocated.'
          call xit("define_regrid",-1)
        endif
        if ( .not. associated(dst_mask_values) ) then
          write(6,*)'define_regrid: Masking dst field but dst_mask_values is not allocated.'
          call xit("define_regrid",-1)
        endif
        write(6,'(2a,i6,a,i6)')"define_regrid: Both src and dst fields are masked. ", &
        "src_mask values = ",src_mask_values,"  dst mask values = ",dst_mask_values
        call ESMF_FieldRegridStore( src_field, dst_field, rc=rc, &
                                   routehandle = regrid_rh_list(id_regrid), &
                                    factorList = factor_list,               &
                               factorIndexList = factor_index_list,         &
                                  srcFracField = src_fracfield,             &
                                  dstFracField = dst_fracfield,             &
                                 srcMaskValues = src_mask_values,           &
                                 dstMaskValues = dst_mask_values,           &
                                  regridmethod = regrid_method,             &
                                    polemethod = pole_method,               &
                               regridPoleNPnts = pole_npts,                 &
                                      lineType = line_type,                 &
                                      normType = norm_type,                 &
                                unmappedaction = unmapped_action,           &
                              ignoreDegenerate = T_esmf,                    &
                             srcTermProcessing = src_term_processing,       &
                                 pipelineDepth = pipeline_depth,            &
                               unmappedDstList = unmapped_dst_list)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Problem creating regrid route handle", line=lnum, method=method)
      else if ( local_use_src_mask ) then
        !--- Mask only the source field
        if ( .not. associated(src_mask_values) ) then
          write(6,*)'define_regrid: Masking src field but src_mask_values is not allocated.'
          call xit("define_regrid",-1)
        endif
        write(6,'(2a,i6)')"define_regrid: Only src field is masked. ", &
                      "src mask values = ",src_mask_values
        call ESMF_FieldRegridStore( src_field, dst_field, rc=rc, &
                                   routehandle = regrid_rh_list(id_regrid), &
                                    factorList = factor_list,               &
                               factorIndexList = factor_index_list,         &
                                  srcFracField = src_fracfield,             &
                                  dstFracField = dst_fracfield,             &
                                 srcMaskValues = src_mask_values,           &
                                  regridmethod = regrid_method,             &
                                    polemethod = pole_method,               &
                               regridPoleNPnts = pole_npts,                 &
                                      lineType = line_type,                 &
                                      normType = norm_type,                 &
                                unmappedaction = unmapped_action,           &
                              ignoreDegenerate = T_esmf,                    &
                             srcTermProcessing = src_term_processing,       &
                                 pipelineDepth = pipeline_depth,            &
                               unmappedDstList = unmapped_dst_list)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Problem creating regrid route handle", line=lnum, method=method)
      else if ( local_use_dst_mask ) then
        !--- Mask only the destination field
        if ( .not. associated(dst_mask_values) ) then
          write(6,*)'define_regrid: Masking dst field but dst_mask_values is not allocated.'
          call xit("define_regrid",-1)
        endif
        write(6,'(2a,i6)')"define_regrid: Only dst field is masked. ", &
                      "dst mask values = ",dst_mask_values
        call ESMF_FieldRegridStore( src_field, dst_field, rc=rc, &
                                   routehandle = regrid_rh_list(id_regrid), &
                                    factorList = factor_list,               &
                               factorIndexList = factor_index_list,         &
                                 dstMaskValues = dst_mask_values,           &
                                  srcFracField = src_fracfield,             &
                                  dstFracField = dst_fracfield,             &
                                  regridmethod = regrid_method,             &
                                    polemethod = pole_method,               &
                               regridPoleNPnts = pole_npts,                 &
                                      lineType = line_type,                 &
                                      normType = norm_type,                 &
                                unmappedaction = unmapped_action,           &
                              ignoreDegenerate = T_esmf,                    &
                             srcTermProcessing = src_term_processing,       &
                                 pipelineDepth = pipeline_depth,            &
                               unmappedDstList = unmapped_dst_list)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Problem creating regrid route handle", line=lnum, method=method)
      else
        !--- No mask used on either field
        write(6,'(a)')"define_regrid: Neither src nor dst fields are masked."
        call ESMF_FieldRegridStore( src_field, dst_field, rc=rc, &
                                   routehandle = regrid_rh_list(id_regrid), &
                                    factorList = factor_list,               &
                               factorIndexList = factor_index_list,         &
                                  srcFracField = src_fracfield,             &
                                  dstFracField = dst_fracfield,             &
                                  regridmethod = regrid_method,             &
                                    polemethod = pole_method,               &
                               regridPoleNPnts = pole_npts,                 &
                                      lineType = line_type,                 &
                                      normType = norm_type,                 &
                                unmappedaction = unmapped_action,           &
                              ignoreDegenerate = T_esmf,                    &
                             srcTermProcessing = src_term_processing,       &
                                 pipelineDepth = pipeline_depth,            &
                               unmappedDstList = unmapped_dst_list)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Problem creating regrid route handle", line=lnum, method=method)
      endif

      !--- Associate the route handle pointer in regrid_list with the
      !--- actual route handle in regrid_rh_list
      regrid_list(id_regrid)%rh => regrid_rh_list(id_regrid)

      !--- Determine the shape of src data array
      if ( associated(farray_ptr) ) nullify( farray_ptr )
      call ESMF_FieldGet( src_field, farrayPtr=farray_ptr, rc=rc )
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Problem getting a pointer to src_field", line=lnum, method=method)
      if ( associated(farray_ptr) ) then
        nx_src = size(farray_ptr,dim=1)
        ny_src = size(farray_ptr,dim=2)
        !--- Do not deallocate any data associated with farray_ptr
        nullify(farray_ptr)
      else
        write(*,*)'Failed to get src field shape for regrid operation.'
        call xit("define_regrid",-1)
      endif

      !--- Determine the shape of dst data array
      if ( associated(farray_ptr) ) nullify( farray_ptr )
      call ESMF_FieldGet( dst_field, farrayPtr=farray_ptr, rc=rc )
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Problem getting a pointer to dst_field", line=lnum, method=method)
      if ( associated(farray_ptr) ) then
        nx_dst = size(farray_ptr,dim=1)
        ny_dst = size(farray_ptr,dim=2)
        !--- Do not deallocate any data associated with farray_ptr
        nullify(farray_ptr)
      else
        write(*,*)'Failed to get dst field shape for regrid operation.'
        call xit("define_regrid",-1)
      endif

      if ( regrid_method == ESMF_REGRIDMETHOD_CONSERVE ) then
        !--- This is a conservative remap
        !--- Create dst_frac for this element of regrid list and copy data from dst_fracfield
        !--- Extract the 2D dst_frac array from dst_fracfield
        if ( associated(farray_ptr) ) nullify( farray_ptr )
        call ESMF_FieldGet( dst_fracfield, farrayPtr=farray_ptr, rc=rc )
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Problem getting a pointer to dst_fracfield", line=lnum, method=method)
        if ( associated(farray_ptr) ) then
          nx_dst = size(farray_ptr,dim=1)
          ny_dst = size(farray_ptr,dim=2)
          if ( associated(regrid_list(id_regrid)%dst_frac) ) then
            deallocate(regrid_list(id_regrid)%dst_frac)
          endif
          allocate( regrid_list(id_regrid)%dst_frac(nx_dst,ny_dst) )
          regrid_list(id_regrid)%dst_frac = farray_ptr
          !--- Do not deallocate any data associated with farray_ptr
          nullify(farray_ptr)
        else
          write(*,*)'Failed to get valid dst fraction pointer from ESMF_FieldGet'
          call xit("define_regrid",-1)
        endif

        !--- Create src_frac for this element of regrid list and copy data from src_fracfield
        !--- Extract the 2D src_frac array from src_fracfield
        if ( associated(farray_ptr) ) nullify( farray_ptr )
        call ESMF_FieldGet( src_fracfield, farrayPtr=farray_ptr, rc=rc )
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Problem getting a pointer to src_fracfield", line=lnum, method=method)
        if ( associated(farray_ptr) ) then
          nx_src = size(farray_ptr,dim=1)
          ny_src = size(farray_ptr,dim=2)
          if ( associated(regrid_list(id_regrid)%src_frac) ) then
            deallocate(regrid_list(id_regrid)%src_frac)
          endif
          allocate( regrid_list(id_regrid)%src_frac(nx_src,ny_src) )
          regrid_list(id_regrid)%src_frac = farray_ptr
          !--- Do not deallocate any data associated with farray_ptr
          nullify(farray_ptr)
        else
          write(*,*)'Failed to get valid src fraction pointer from ESMF_FieldGet'
          call xit("define_regrid",-1)
        endif

      else
        !--- This is not a conservative remap

        !--- In this case assign src_frac and dst_frac with src and dst mask values
        !--- found in the associated ESMF grid objects
        if ( associated(cpl_grid(id_src_grid)%mask) ) then
          if ( nx_src /= size(cpl_grid(id_src_grid)%mask,dim=1) .or. &
               ny_src /= size(cpl_grid(id_src_grid)%mask,dim=2) ) then
            write(6,*)"define_regrid: Shape mismatch between SRC field and grid mask"
            write(6,*)"  shape(field) = ",nx_src,ny_src
            write(6,*)"  shape( mask) = ",shape(cpl_grid(id_src_grid)%mask)
            call xit("define_regrid",-1)
          endif
          if ( associated(regrid_list(id_regrid)%src_frac) ) then
            deallocate(regrid_list(id_regrid)%src_frac)
          endif
          allocate( regrid_list(id_regrid)%src_frac(nx_src,ny_src) )
          regrid_list(id_regrid)%src_frac = 1
          if ( associated(src_mask_values) ) then
            do idx1=1,size(src_mask_values)
              where (cpl_grid(id_src_grid)%mask == src_mask_values(idx1)) &
                    regrid_list(id_regrid)%src_frac = 0
            enddo
          endif
        else
          write(*,*)'define_regrid: mask is not associated with SRC grid ',id_src_grid
          call xit("define_regrid",-1)
        endif

        if ( associated(cpl_grid(id_dst_grid)%mask) ) then
          if ( nx_dst /= size(cpl_grid(id_dst_grid)%mask,dim=1) .or. &
               ny_dst /= size(cpl_grid(id_dst_grid)%mask,dim=2) ) then
            write(6,*)"define_regrid: Shape mismatch between DST field and grid mask"
            write(6,*)"  shape(field) = ",nx_dst,ny_dst
            write(6,*)"  shape( mask) = ",shape(cpl_grid(id_dst_grid)%mask)
            call xit("define_regrid",-1)
          endif
          if ( associated(regrid_list(id_regrid)%dst_frac) ) then
            deallocate(regrid_list(id_regrid)%dst_frac)
          endif
          allocate( regrid_list(id_regrid)%dst_frac(nx_dst,ny_dst) )
          regrid_list(id_regrid)%dst_frac = 1
          if ( associated(dst_mask_values) ) then
            do idx1=1,size(dst_mask_values)
              where (cpl_grid(id_dst_grid)%mask == dst_mask_values(idx1)) &
                    regrid_list(id_regrid)%dst_frac = 0
            enddo
          endif
        else
          write(*,*)'define_regrid: mask is not associated with DST grid ',id_dst_grid
          call xit("define_regrid",-1)
        endif

      endif

      !--- In all cases extract the src and dst grid areas from the src and dst ESMF grid objects

      !--- Retrieve a pointer to the grid area in the src grid object
      nullify(src_area)
      call ESMF_GridGetItem(src_grid, &
                            staggerloc=ESMF_STAGGERLOC_CENTER, &
                            itemflag=ESMF_GRIDITEM_AREA, &
                            farrayPtr=src_area, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot retrieve src grid area", line=lnum, method=method)

      !--- Assign grid cell areas with areas from the associated ESMF grid objects
      if ( associated(src_area) ) then
        if ( nx_src /= size(src_area,dim=1) .or. &
             ny_src /= size(src_area,dim=2) ) then
          write(6,*)"define_regrid: Shape mismatch between SRC field and grid area"
          write(6,*)"  shape(field) = ",nx_src,ny_src
          write(6,*)"  shape( area) = ",shape(src_area)
          call xit("define_regrid",-1)
        endif
        if ( associated( regrid_list(id_regrid)%src_area ) ) then
          deallocate(regrid_list(id_regrid)%src_area)
        endif
        allocate( regrid_list(id_regrid)%src_area(nx_src,ny_src) )
        regrid_list(id_regrid)%src_area = src_area
        !--- Do not deallocate src_area
        nullify(src_area)
      else
        write(*,*)'define_regrid: Area is not associated with SRC grid ',id_src_grid
        call xit("define_regrid",-1)
      endif

      !--- Retrieve a pointer to the grid area in the dst grid object
      nullify(dst_area)
      call ESMF_GridGetItem(dst_grid, &
                            staggerloc=ESMF_STAGGERLOC_CENTER, &
                            itemflag=ESMF_GRIDITEM_AREA, &
                            farrayPtr=dst_area, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot retrieve dst grid area", line=lnum, method=method)

      if ( associated(dst_area) ) then
        if ( nx_dst /= size(dst_area,dim=1) .or. &
             ny_dst /= size(dst_area,dim=2) ) then
          write(6,*)"define_regrid: Shape mismatch between DST field and grid area"
          write(6,*)"  shape(field) = ",nx_dst,ny_dst
          write(6,*)"  shape( area) = ",shape(dst_area)
          call xit("define_regrid",-1)
        endif
        if ( associated( regrid_list(id_regrid)%dst_area ) ) then
          deallocate(regrid_list(id_regrid)%dst_area)
        endif
        allocate( regrid_list(id_regrid)%dst_area(nx_dst,ny_dst) )
        regrid_list(id_regrid)%dst_area = dst_area
        !--- Do not deallocate dst_area
        nullify(dst_area)
      else
        write(*,*)'define_regrid: Area is not associated with DST grid ',id_dst_grid
        call xit("define_regrid",-1)
      endif

      call cpu_time( time_out )
      time_out = time_out - time_in
      write(6,'(a,g15.7)')'define_regrid: Time in ESMF_FieldRegridStore is ',time_out
      call flush(6)

      if ( associated(unmapped_dst_list) ) then
        !--- There are unmapped destination points
        n_unmapped = size(unmapped_dst_list)

        if ( verbose > 0 ) then
          write(6,'(a,i6,a)')'define_regrid: There are ', &
              n_unmapped,' unmapped destination grid cells.'
          call flush(6)
        endif

        if ( (verbose > 1 .or. n_unmapped < 31) .and. n_unmapped > 0 ) then
          !--- Print the list of unmapped destination points

          !--- Get an ESMF array containing dim 1 (lon) center coordinate values
          call ESMF_GridGetCoord(dst_grid, coordDim=one, array=coordAry, &
                                 staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc)
          lnum = __LINE__ ; call handle_ESMF_err(rc, &
              "Unable to retrieve dst_grid coord array", line=lnum, method=method)

          !--- Determine the rank of the ESMF array
          call ESMF_ArrayGet(coordAry, rank=rank, rc=rc)
          lnum = __LINE__ ; call handle_ESMF_err(rc, &
              "Unable to retrieve rank of dst_grid coord array", line=lnum, method=method)

          !--- Determine dimension sizes from dst_grid according to rank
          if ( rank == 1 ) then
            !--- 1st coord is 1D and depends on 1st Grid dim
            !--- 2nd coord is 1D and depends on 2nd Grid dim
            call ESMF_GridGetCoordBounds(dst_grid, coordDim=one,  &
                                         exclusiveCount=excCount, &
                                         staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc)
            nx_c = excCount(1)
            call ESMF_GridGetCoordBounds(dst_grid, coordDim=two,  &
                                         exclusiveCount=excCount, &
                                         staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc)
            ny_c = excCount(1)
            call ESMF_GridGetCoordBounds(dst_grid, coordDim=one,  &
                                         exclusiveCount=excCount, &
                                         staggerloc=ESMF_STAGGERLOC_CORNER, rc=rc)
            nx_v = excCount(1)
            call ESMF_GridGetCoordBounds(dst_grid, coordDim=two,  &
                                         exclusiveCount=excCount, &
                                         staggerloc=ESMF_STAGGERLOC_CORNER, rc=rc)
            ny_v = excCount(1)
          else if ( rank == 2 ) then
            !--- 1st coord is 2D and depends on 1st and 2nd Grid dim
            !--- 2nd coord is 2D and depends on 1st and 2nd Grid dim
            !--- Assume both dimensions have the same shape
            call ESMF_GridGetCoordBounds(dst_grid, coordDim=one,  &
                                         exclusiveCount=excCount, &
                                         staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc)
            nx_c = excCount(1)
            ny_c = excCount(2)
            call ESMF_GridGetCoordBounds(dst_grid, coordDim=one,  &
                                         exclusiveCount=excCount, &
                                         staggerloc=ESMF_STAGGERLOC_CORNER, rc=rc)
            nx_v = excCount(1)
            ny_v = excCount(2)
          else
            write(*,*)'print_grid_info: Invalid rank for dst_grid coord array = ',rank
            call xit("define_regrid",-1)
          endif
          if ( verbose > 10 ) then
            write(6,*)"define_regrid: nx_c, ny_c = ",nx_c,ny_c
            write(6,*)"define_regrid: nx_v, ny_v = ",nx_v,ny_v
          endif

          !--- Extract lon_c coordinate values from this array
          if ( rank == 1 ) then
            call ESMF_ArrayGet(coordAry, farrayPtr=coord1d, rc=rc)
            lnum = __LINE__ ; call handle_ESMF_err(rc, &
                "Unable to retrieve 1d pointer to dst_grid coord array", &
                line=lnum, method=method)
            if ( size(coord1d) /= nx_c ) then
              write(6,*)"define_regrid: Size mismatch for lon_c"
              write(6,*)"    found nx = ",shape(coord1d)
              write(6,*)" expected nx = ",nx_c
              call xit("define_regrid",-1)
            endif
            allocate( lon_c(nx_c,ny_c) )
            do idx1=1,ny_c
              lon_c(1:nx_c,idx1) = coord1d(1:nx_c)
            enddo
            nullify(coord1d)
          else if ( rank == 2 ) then
            call ESMF_ArrayGet(coordAry, farrayPtr=coord2d, rc=rc)
            lnum = __LINE__ ; call handle_ESMF_err(rc, &
                "Unable to retrieve 2d pointer to dst_grid coord array", &
                line=lnum, method=method)
            if ( size(coord2d,dim=1) /= nx_c .or. size(coord2d,dim=2) /= ny_c ) then
              write(6,*)"define_regrid: Size mismatch for lon_c"
              write(6,*)"    found nx, ny = ",shape(coord2d)
              write(6,*)" expected nx, ny = ",nx_c,ny_c
              call xit("define_regrid",-1)
            endif
            allocate( lon_c(nx_c,ny_c) )
            lon_c = coord2d
            nullify(coord2d)
          else
            write(*,*)'print_grid_info: Invalid rank for dst_grid coord array = ',rank
            call xit("define_regrid",-1)
          endif

          !--- Get an ESMF array containing dim 2 (lat) center coordinate values
          call ESMF_GridGetCoord(dst_grid, coordDim=two, array=coordAry, &
                                 staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc)
          lnum = __LINE__ ; call handle_ESMF_err(rc, &
              "Unable to retrieve dst_grid coord array", line=lnum, method=method)

          !--- Extract lat_c coordinate values from this array
          call ESMF_ArrayGet(coordAry, rank=rank, rc=rc)
          lnum = __LINE__ ; call handle_ESMF_err(rc, &
              "Unable to retrieve rank of dst_grid coord array", line=lnum, method=method)
          if ( rank == 1 ) then
            call ESMF_ArrayGet(coordAry, farrayPtr=coord1d, rc=rc)
            lnum = __LINE__ ; call handle_ESMF_err(rc, &
                "Unable to retrieve 1d pointer to dst_grid coord array", &
                line=lnum, method=method)
            if ( size(coord1d) /= ny_c ) then
              write(6,*)"define_regrid: Size mismatch for lat_c"
              write(6,*)"    found ny = ",shape(coord1d)
              write(6,*)" expected ny = ",ny_c
              call xit("define_regrid",-1)
            endif
            allocate( lat_c(nx_c,ny_c) )
            do idx1=1,nx_c
              lat_c(idx1,1:ny_c) = coord1d(1:ny_c)
            enddo
            nullify(coord1d)
          else if ( rank == 2 ) then
            call ESMF_ArrayGet(coordAry, farrayPtr=coord2d, rc=rc)
            lnum = __LINE__ ; call handle_ESMF_err(rc, &
                "Unable to retrieve 2d pointer to dst_grid coord array", &
                line=lnum, method=method)
            if ( size(coord2d,dim=1) /= nx_c .or. size(coord2d,dim=2) /= ny_c ) then
              write(6,*)"define_regrid: Size mismatch for lat_c"
              write(6,*)"    found nx, ny = ",shape(coord2d)
              write(6,*)" expected nx, ny = ",nx_c,ny_c
              call xit("define_regrid",-1)
            endif
            allocate( lat_c(nx_c,ny_c) )
            lat_c = coord2d
            nullify(coord2d)
          else
            write(*,*)'print_grid_info: Invalid rank for dst_grid coord array = ',rank
            call xit("define_regrid",-1)
          endif

          !--- Get an ESMF array containing dim 1 (lon) corner coordinate values
          call ESMF_GridGetCoord(dst_grid, coordDim=one, array=coordAry, &
                                 staggerloc=ESMF_STAGGERLOC_CORNER, rc=rc)
          lnum = __LINE__ ; call handle_ESMF_err(rc, &
              "Unable to retrieve dst_grid coord array", line=lnum, method=method)

          !--- Extract lon_v coordinate values from this array
          call ESMF_ArrayGet(coordAry, rank=rank, rc=rc)
          lnum = __LINE__ ; call handle_ESMF_err(rc, &
              "Unable to retrieve rank of dst_grid coord array", line=lnum, method=method)
          if ( rank == 1 ) then
            call ESMF_ArrayGet(coordAry, farrayPtr=coord1d, rc=rc)
            lnum = __LINE__ ; call handle_ESMF_err(rc, &
                "Unable to retrieve 1d pointer to dst_grid coord array", &
                line=lnum, method=method)
            if ( size(coord1d) /= nx_v ) then
              write(6,*)"define_regrid: Size mismatch for lon_v"
              write(6,*)"    found ny = ",shape(coord1d)
              write(6,*)" expected ny = ",nx_v
              call xit("define_regrid",-1)
            endif
            allocate( lon_v(nx_v,ny_v) )
            do idx1=1,ny_v
              lon_v(1:nx_v,idx1) = coord1d(1:nx_v)
            enddo
            nullify(coord1d)
          else if ( rank == 2 ) then
            call ESMF_ArrayGet(coordAry, farrayPtr=coord2d, rc=rc)
            lnum = __LINE__ ; call handle_ESMF_err(rc, &
                "Unable to retrieve 2d pointer to dst_grid coord array", &
                line=lnum, method=method)
            if ( size(coord2d,dim=1) /= nx_v .or. size(coord2d,dim=2) /= ny_v ) then
              write(6,*)"define_regrid: Size mismatch for lon_v"
              write(6,*)"    found nx, ny = ",shape(coord2d)
              write(6,*)" expected nx, ny = ",nx_v,ny_v
              call xit("define_regrid",-1)
            endif
            allocate( lon_v(nx_v,ny_v) )
            lon_v = coord2d
            nullify(coord2d)
          else
            write(*,*)'print_grid_info: Invalid rank for dst_grid coord array = ',rank
            call xit("define_regrid",-1)
          endif

          !--- Get an ESMF array containing dim 2 (lat) corner coordinate values
          call ESMF_GridGetCoord(dst_grid, coordDim=two, array=coordAry, &
                                 staggerloc=ESMF_STAGGERLOC_CORNER, rc=rc)
          lnum = __LINE__ ; call handle_ESMF_err(rc, &
              "Unable to retrieve dst_grid coord array", line=lnum, method=method)

          !--- Extract lat_v coordinate values from this array
          call ESMF_ArrayGet(coordAry, rank=rank, rc=rc)
          lnum = __LINE__ ; call handle_ESMF_err(rc, &
              "Unable to retrieve rank of dst_grid coord array", line=lnum, method=method)
          if ( rank == 1 ) then
            call ESMF_ArrayGet(coordAry, farrayPtr=coord1d, rc=rc)
            lnum = __LINE__ ; call handle_ESMF_err(rc, &
                "Unable to retrieve 1d pointer to dst_grid coord array", &
                line=lnum, method=method)
            if ( size(coord1d) /= ny_v ) then
              write(6,*)"define_regrid: Size mismatch for lat_v"
              write(6,*)"    found ny = ",shape(coord1d)
              write(6,*)" expected ny = ",ny_v
              call xit("define_regrid",-1)
            endif
            allocate( lat_v(nx_v,ny_v) )
            do idx1=1,nx_v
              lat_v(idx1,1:ny_v) = coord1d(1:ny_v)
            enddo
            nullify(coord1d)
          else if ( rank == 2 ) then
            call ESMF_ArrayGet(coordAry, farrayPtr=coord2d, rc=rc)
            lnum = __LINE__ ; call handle_ESMF_err(rc, &
                "Unable to retrieve 2d pointer to dst_grid coord array", &
                line=lnum, method=method)
            if ( size(coord2d,dim=1) /= nx_v .or. size(coord2d,dim=2) /= ny_v ) then
              write(6,*)"define_regrid: Size mismatch for lat_v"
              write(6,*)"    found nx, ny = ",shape(coord2d)
              write(6,*)" expected nx, ny = ",nx_v,ny_v
              call xit("define_regrid",-1)
            endif
            allocate( lat_v(nx_v,ny_v) )
            lat_v = coord2d
            nullify(coord2d)
          else
            write(*,*)'print_grid_info: Invalid rank for dst_grid coord array = ',rank
            call xit("define_regrid",-1)
          endif

          do idx1=1,n_unmapped
            !--- Determine indicies on the 2D grid assuming column major order
            idx2 = unmapped_dst_list(idx1)
            ix = 1 + mod(idx2-1,nx_c)
            iy = (idx2-1)/nx_c + 1
            ixp1 = ix+1
            iyp1 = iy+1
            !--- Assume cyclic longitude when ix+1 > nx_v
            if ( ixp1 > nx_v ) ixp1 = ixp1 - nx_v
            write(6,'(i6,"  index = ",i6,"  lon, lat = ",2f15.10,$)') &
               idx1, idx2, lon_c(ix,iy), lat_c(ix,iy)
            !--- The corners are at (ix,iy), (ixp1,iy), (ixp1,iyp1), (ix,iyp1)
            write(6,'(a,4(" (",f12.7,",",f12.7,")"))') &
               "  corners:", &
               lon_v(ix,iy),     lat_v(ix,iy),     &
               lon_v(ixp1,iy),   lat_v(ixp1,iy),   &
               lon_v(ixp1,iyp1), lat_v(ixp1,iyp1), &
               lon_v(ix,iyp1),   lat_v(ix,iyp1)
          enddo
          call flush(6)

          deallocate( lat_c, lon_c, lat_v, lon_v )

          if ( .not. local_ignore_unmapped) then
            !--- Abort if there are unmapped points
            call xit("UNMAPPED",-1)
          endif

        endif
      else
        n_unmapped = 0
      endif

      !--- Reformat source and destination grid area and grid fraction as 1D arrays
      !--- that are then placed in internal storage. These same arrays are written
      !--- to the weights file if it is created below.
      if ( regrid_method == ESMF_REGRIDMETHOD_CONSERVE ) then
        !--- This is a conservative remap

        !--- Compute src grid cell areas and store them in the src_field object
        call ESMF_FieldRegridGetArea(src_field, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Problem getting src_field area", line=lnum, method=method)

        !--- Extract the 2D area array from src_field and assign it to a 1D array
        if ( associated(farray_ptr) ) nullify( farray_ptr )
        call ESMF_FieldGet( src_field, farrayPtr=farray_ptr, rc=rc )
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Problem getting src_field area pointer", line=lnum, method=method)
        if ( associated(farray_ptr) ) then
          nx = size(farray_ptr,dim=1)
          ny = size(farray_ptr,dim=2)
          if ( verbose > 10 ) then
            write(6,'(a,2i8)')'define_regrid: conserve src_field area nx,ny = ',nx,ny
            call flush(6)
          endif

          allocate( srcArea_out(nx*ny) )
          srcArea_out = reshape( farray_ptr, (/nx*ny/) )
          nullify( farray_ptr )

          !--- Add srcArea_out to internal storage
          strng = " "
          write(strng,'(a,i3.3)')"src_area",id_regrid

          !--- Define a cpl_vinfo_t variable to use for src_area
          allocate( tmp_vinfo )
          tmp_vinfo%name     = trim(strng)
          tmp_vinfo%grid     = "area_grid"
          tmp_vinfo%tag      = cpl_var_new_tag(5000)
          tmp_vinfo%size     = nx*ny
          tmp_vinfo%nlon     = nx  !--- Do not include overlap
          tmp_vinfo%nlat     = ny
          tmp_vinfo%olap     = 0
          tmp_vinfo%ncopy    = 1
          tmp_vinfo%freq     = 0
          tmp_vinfo%ccc_name = "AREA"
          tmp_vinfo%regrid%method   = trim(map_method)
          tmp_vinfo%regrid%norm_opt = trim(this_norm_type)
          tmp_vinfo%regrid%masked   = local_use_src_mask
          !--- Add this variable defintion
          call add_cpl_var(tmp_vinfo, exists=0)
          nullify( tmp_vinfo )

          !--- The src_area record will be created if it does not exist
          call store_AppendRecord(trim(strng), other_gid, srcArea_out, create=.true.)
          call store_ListRecord(trim(strng), other_gid)

        else
          write(*,*)'Failed to get valid cell area pointer from ESMF_FieldGet'
          call xit("define_regrid",-1)
        endif

        !--- Compute dst grid cell areas and store them in the dst_field object
        call ESMF_FieldRegridGetArea(dst_field, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Problem getting dst_field area", line=lnum, method=method)

        !--- Extract the 2D area array from dst_field and assign it to a 1D array
        if ( associated(farray_ptr) ) nullify( farray_ptr )
        call ESMF_FieldGet( dst_field, farrayPtr=farray_ptr, rc=rc )
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Problem getting dst_field area pointer", line=lnum, method=method)
        if ( associated(farray_ptr) ) then
          nx = size(farray_ptr,dim=1)
          ny = size(farray_ptr,dim=2)
          if ( verbose > 10 ) then
            write(6,'(a,2i8)')'define_regrid: conserve dst_field area nx,ny = ',nx,ny
            call flush(6)
          endif

          allocate( dstArea_out(nx*ny) )
          dstArea_out = reshape( farray_ptr, (/nx*ny/) )
          nullify( farray_ptr )

          !--- Add dstArea_out to internal storage
          strng = " "
          write(strng,'(a,i3.3)')"dst_area",id_regrid

          !--- Define a cpl_vinfo_t variable to use for dst_area
          allocate( tmp_vinfo )
          tmp_vinfo%name     = trim(strng)
          tmp_vinfo%grid     = "area_grid"
          tmp_vinfo%tag      = cpl_var_new_tag(5000)
          tmp_vinfo%size     = nx*ny
          tmp_vinfo%nlon     = nx  !--- Do not include overlap
          tmp_vinfo%nlat     = ny
          tmp_vinfo%olap     = 0
          tmp_vinfo%ncopy    = 1
          tmp_vinfo%freq     = 0
          tmp_vinfo%ccc_name = "AREA"
          tmp_vinfo%regrid%method   = trim(map_method)
          tmp_vinfo%regrid%norm_opt = trim(this_norm_type)
          tmp_vinfo%regrid%masked   = local_use_dst_mask
          !--- Add this variable defintion
          call add_cpl_var(tmp_vinfo, exists=0)
          nullify( tmp_vinfo )

          !--- The dst_area record will be created if it does not exist
          call store_AppendRecord(trim(strng), other_gid, dstArea_out, create=.true.)
          call store_ListRecord(trim(strng), other_gid)

        else
          write(*,*)'Failed to get valid cell area pointer from ESMF_FieldGet'
          call xit("define_regrid",-1)
        endif

        !--- Copy the 2D src_fracfield array to the 1D srcFrac_out array
        if ( associated(farray_ptr) ) nullify( farray_ptr )
        call ESMF_FieldGet( src_fracfield, farrayPtr=farray_ptr, rc=rc )
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Problem getting src_fracfield", line=lnum, method=method)
        if ( associated(farray_ptr) ) then
          nx = size(farray_ptr,dim=1)
          ny = size(farray_ptr,dim=2)
          if ( verbose > 10 ) then
            write(6,'(a,2i8)')'define_regrid: conserve src_fracfield nx,ny = ',nx,ny
            call flush(6)
          endif
          allocate( srcFrac_out(nx*ny) )
          srcFrac_out = reshape( farray_ptr, (/nx*ny/) )
          nullify( farray_ptr )

          !--- Add srcFrac_out to internal storage
          strng = " "
          write(strng,'(a,i3.3)')"src_frac",id_regrid

          !--- Define a cpl_vinfo_t variable to use for src_frac
          allocate( tmp_vinfo )
          tmp_vinfo%name     = trim(strng)
          tmp_vinfo%grid     = "frac_grid"
          tmp_vinfo%tag      = cpl_var_new_tag(5000)
          tmp_vinfo%size     = nx*ny
          tmp_vinfo%nlon     = nx  !--- Do not include overlap
          tmp_vinfo%nlat     = ny
          tmp_vinfo%olap     = 0
          tmp_vinfo%ncopy    = 1
          tmp_vinfo%freq     = 0
          tmp_vinfo%ccc_name = "SFRC"
          tmp_vinfo%regrid%method   = trim(map_method)
          tmp_vinfo%regrid%norm_opt = trim(this_norm_type)
          tmp_vinfo%regrid%masked   = local_use_src_mask
          !--- Add this variable defintion
          call add_cpl_var(tmp_vinfo, exists=0)
          nullify( tmp_vinfo )

          !--- The src_frac record will be created if it does not exist
          call store_AppendRecord(trim(strng), other_gid, srcFrac_out, create=.true.)
          call store_ListRecord(trim(strng), other_gid)

        else
          write(*,*)'Failed to get valid src_fracfield pointer from ESMF_FieldGet'
          call xit("define_regrid",-1)
        endif

        !--- Copy the 2D dst_fracfield array to the 1D dstFrac_out array
        if ( associated(farray_ptr) ) nullify( farray_ptr )
        call ESMF_FieldGet( dst_fracfield, farrayPtr=farray_ptr, rc=rc )
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Problem getting dst_fracfield", line=lnum, method=method)
        if ( associated(farray_ptr) ) then
          nx = size(farray_ptr,dim=1)
          ny = size(farray_ptr,dim=2)
          if ( verbose > 10 ) then
            write(6,'(a,2i8)')'define_regrid: other dst_fracfield nx,ny = ',nx,ny
            call flush(6)
          endif
          allocate( dstFrac_out(nx*ny) )
          dstFrac_out = reshape( farray_ptr, (/nx*ny/) )
          nullify( farray_ptr )

          !--- Add dstFrac_out to internal storage
          strng = " "
          write(strng,'(a,i3.3)')"dst_frac",id_regrid

          !--- Define a cpl_vinfo_t variable to use for dst_frac
          allocate( tmp_vinfo )
          tmp_vinfo%name     = trim(strng)
          tmp_vinfo%grid     = "frac_grid"
          tmp_vinfo%tag      = cpl_var_new_tag(5000)
          tmp_vinfo%size     = nx*ny
          tmp_vinfo%nlon     = nx  !--- Do not include overlap
          tmp_vinfo%nlat     = ny
          tmp_vinfo%olap     = 0
          tmp_vinfo%ncopy    = 1
          tmp_vinfo%freq     = 0
          tmp_vinfo%ccc_name = "DFRC"
          tmp_vinfo%regrid%method   = trim(map_method)
          tmp_vinfo%regrid%norm_opt = trim(this_norm_type)
          tmp_vinfo%regrid%masked   = local_use_dst_mask
          !--- Add this variable defintion
          call add_cpl_var(tmp_vinfo, exists=0)
          nullify( tmp_vinfo )

          !--- The dst_frac record will be created if it does not exist
          call store_AppendRecord(trim(strng), other_gid, dstFrac_out, create=.true.)
          call store_ListRecord(trim(strng), other_gid)
        else
          write(*,*)'Failed to get valid dst_fracfield pointer from ESMF_FieldGet'
          call xit("define_regrid",-1)
        endif

      endif

#ifdef save_regrid_weights
      !--- Optionally save weights to an external file
      save_weights_to_file = .true.
      if ( save_weights_to_file ) then

        !--- Save these weights to an external file
        write(wgt_file,'(a,i3.3,a)')"regrid_weights",id_regrid,".nc"
        regrid_list(id_regrid)%wgt_grid_file_name = trim(wgt_file)
        if (verbose > 10) then
          write(6,'(2a)')"define_regrid: Creating weights file ",trim(wgt_file)
          call flush(6)
        endif

        grid_file_type = grid_desc_file_type(src_grid_file_name)
        select case (grid_file_type)
          case (1) !--- A SCRIP format file
            call read_scrip_grid_file(src_grid_file_name, title=src_grid_title)
          case (2) !---A NEMO "mesh_mask" format file
            call read_NEMO_mesh_mask_file(src_grid_file_name, title=src_grid_title)
          case default
            write(6,*)'define_regrid: ',trim(src_grid_file_name), &
                      '  Invalid grid_file_type = ',grid_file_type
            call xit("define_regrid",-1)
        end select

        grid_file_type = grid_desc_file_type(dst_grid_file_name)
        select case (grid_file_type)
          case (1) !--- A SCRIP format file
            call read_scrip_grid_file(dst_grid_file_name, title=dst_grid_title)
          case (2) !---A NEMO "mesh_mask" format file
            call read_NEMO_mesh_mask_file(dst_grid_file_name, title=dst_grid_title)
          case default
            write(6,*)'define_regrid: ',trim(dst_grid_file_name), &
                      '  Invalid grid_file_type = ',grid_file_type
            call xit("define_regrid",-1)
        end select

        wgt_title = " "
        if ( regrid_method == ESMF_REGRIDMETHOD_CONSERVE ) then
          write(wgt_title,'(4a)')"Weights for conservative remap from ", &
              trim(src_grid_title)," to ",trim(dst_grid_title)
        else if ( regrid_method == ESMF_REGRIDMETHOD_BILINEAR ) then
          write(wgt_title,'(4a)')"Weights for bilinear remap from ", &
              trim(src_grid_title)," to ",trim(dst_grid_title)
        else if ( regrid_method == ESMF_REGRIDMETHOD_PATCH ) then
          write(wgt_title,'(4a)')"Weights for patch recovery remap from ", &
              trim(src_grid_title)," to ",trim(dst_grid_title)
        else if ( regrid_method == ESMF_REGRIDMETHOD_NEAREST_STOD ) then
          write(wgt_title,'(4a)')"Weights for nearest src to dst remap from ", &
              trim(src_grid_title)," to ",trim(dst_grid_title)
        else if ( regrid_method == ESMF_REGRIDMETHOD_NEAREST_DTOS ) then
          write(wgt_title,'(4a)')"Weights for nearest dst to src remap from ", &
              trim(src_grid_title)," to ",trim(dst_grid_title)
        else
          !--- Do not consider this an error
          write(wgt_title,'(4a)')"Weights for remap from ", &
              trim(src_grid_title)," to ",trim(dst_grid_title)
        endif
        if (verbose > 0) then
          write(6,'(2a)')"define_regrid: Weights title: ",trim(wgt_title)
          call flush(6)
        endif

        if ( regrid_method == ESMF_REGRIDMETHOD_CONSERVE ) then
          !--- This is a conservative remap

          !--- Create the external weights file
          call ESMF_OutputScripWeightFile (wgt_file, factor_List, factor_Index_List, &
                             srcFile = trim(src_grid_file_name), &
                             dstFile = trim(dst_grid_file_name), &
                         srcFileType = ESMF_FILEFORMAT_SCRIP,    &
                         dstFileType = ESMF_FILEFORMAT_SCRIP,    &
                              method = regrid_method,            &
                            normType = norm_type,                &
                               title = trim(wgt_title),          &
                             srcArea = srcArea_out,              &
                             dstArea = dstArea_out,              &
                             srcFrac = srcFrac_out,              &
                             dstFrac = dstFrac_out,              &
                                  rc = rc)
          lnum = __LINE__ ; call handle_ESMF_err(rc, &
              "Unable to write weights file", line=lnum, method=method)

        else
          !--- This is a non-conservative remap

          !--- Create the external weights file
          call ESMF_OutputScripWeightFile (wgt_file, factor_List, factor_Index_List, &
                             srcFile = trim(src_grid_file_name), &
                             dstFile = trim(dst_grid_file_name), &
                         srcFileType = ESMF_FILEFORMAT_SCRIP,    &
                         dstFileType = ESMF_FILEFORMAT_SCRIP,    &
                              method = regrid_method,            &
                            normType = norm_type,                &
                               title = trim(wgt_title),          &
                                  rc = rc)
          lnum = __LINE__ ; call handle_ESMF_err(rc, &
              "Unable to write weights file", line=lnum, method=method)

        endif

!xxx ESMF_RegridWeightGenCheck requires the following modules
!xxx use ESMF_RegridWeightGenMod
!xxx use ESMF_RegridWeightGenCheckMod
!xxx          !--- Print some info generated from this weights file by interpolating
!xxx          !--- a known function namely (2 + cos(lat)**2 * cos(2*lon))
!xxx          call ESMF_RegridWeightGenCheck(wgt_file, rc=rc)
!xxx          lnum = __LINE__ ; call handle_ESMF_err(rc, &
!xxx              "Problem checking weights", line=lnum, method=method)

        if ( verbose > 0 ) then
          write(6,'(2a)')'define_regrid: Created weights file ',trim(wgt_file)
          call flush(6)
        endif

      endif
#endif

      !--- Release resources used by this routine

      !--- Mask values used when calculating weights
      if ( associated(src_mask_values) ) deallocate(src_mask_values)
      if ( associated(dst_mask_values) ) deallocate(dst_mask_values)

      !--- These may be allocated when source and destination grid area
      !--- and grid fraction are written to internal storage
      if ( associated(srcFrac_out) ) deallocate( srcFrac_out )
      if ( associated(dstFrac_out) ) deallocate( dstFrac_out )

      !--- unmapped_dst_list may be allocated during the call to ESMF_FieldRegridStore
      if ( associated(unmapped_dst_list) ) deallocate( unmapped_dst_list )

      !--- These are allocated during the call to ESMF_FieldRegridStore
      if ( associated(factor_list) )       deallocate( factor_list )
      if ( associated(factor_index_list) ) deallocate( factor_index_list )

      !--- Release resources used by src_fracfield and dst_fracfield
      call ESMF_FieldDestroy(src_fracfield, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to destroy src_fracfield", line=lnum, method=method)

      call ESMF_FieldDestroy(dst_fracfield, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to destroy dst_fracfield", line=lnum, method=method)

      !--- Release resources used by src_field and dst_field
      call ESMF_FieldDestroy(src_field, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to destroy src_field", line=lnum, method=method)

      call ESMF_FieldDestroy(dst_field, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to destroy dst_field", line=lnum, method=method)

    end subroutine define_regrid

    !***************************************************************************
    !--- Wrapper for a conservative, fracarea remap
    !***************************************************************************
    subroutine regrid_field_cv(src_name, src_data, dst_name, dst_data, alt_name)
      character(*), intent(in) :: src_name, dst_name
      real(kind=8), intent(in) :: src_data(:)
      real(kind=8), intent(out), pointer :: dst_data(:)
      character(*), intent(in), optional :: alt_name

      !--- Local
      character(128) :: src_grid, dst_grid
      integer :: status
      character(64) :: src_fld_name
      type(cpl_vinfo_t), pointer :: curr_vinfo => null()

      if ( present(alt_name) ) then
        src_fld_name = trim(alt_name)
      else
        src_fld_name = trim(src_name)
      endif

      nullify(curr_vinfo)
      allocate(curr_vinfo)

      src_grid = " "
      curr_vinfo = find_cpl_vinfo( name=trim(src_name), status=status )
      if ( status == 0 ) then
        !--- The variable was found in cpl_var_list
        src_grid = trim(curr_vinfo%grid)
      else
        !--- This variable is unknown (ie is not in cpl_var_list)
        write(6,*)"regrid_field_cv: Missing variable ",trim(src_name)
        call xit("regrid_field_cv",-1)
      endif

      dst_grid = " "
      curr_vinfo = find_cpl_vinfo( name=trim(dst_name), status=status )
      if ( status == 0 ) then
        !--- The variable was found in cpl_var_list
        dst_grid = trim(curr_vinfo%grid)
      else
        !--- This variable is unknown (ie is not in cpl_var_list)
        write(6,*)"regrid_field_cv: Missing variable ",trim(dst_name)
        call xit("regrid_field_cv",-2)
      endif

      !--- Always mask the NEMO land points, whether interpolating atm->ocn or ocn->atm
      !--- Cells that do not take part in the remap are assigned a zero value
      if ( trim(adjustl(src_grid)) .eq. "atm" .and. &
           trim(adjustl(dst_grid)) .eq. "ocn") then 
        !--- Conservative-fracarea remap from AGCM to NEMO grid with land values on
        !--- the NEMO grid masked (land does not take part in the remap calculation)
        call regrid_Field(src_name, src_data, dst_name, dst_data, &
                          method="conservative", norm_type="fracarea", &
                          src_fld_name=trim(src_fld_name), &
                          dst_mask_values=(/0/), missing_value=0.0_8)

      else if ( trim(adjustl(src_grid)) .eq. "ocn" .and. &
                trim(adjustl(dst_grid)) .eq. "atm") then 
        !--- Conservative-fracarea remap from NEMO to AGCM grid with land values on
        !--- the NEMO grid masked (land does not take part in the remap calculation)
        call regrid_Field(src_name, src_data, dst_name, dst_data, &
                          method="conservative", norm_type="fracarea", &
                          src_fld_name=trim(src_fld_name), &
                          src_mask_values=(/0/), missing_value=0.0_8)

      else
        !--- The remap is not yet supported
        write(6,*)"regrid_field_cv: Remap from ",trim(src_grid)," grid to ", &
                  trim(dst_grid)," is not currently supported by regrid_field_cv"
        call xit("regrid_field_cv",-3)
      endif

      !--- Clean up
      deallocate(curr_vinfo)

    end subroutine regrid_field_cv

    !***************************************************************************
    !--- Regrid one field onto another field
    !--- NOTE: dst_field will always be reallocated to be the size of dst_name
    !***************************************************************************
    subroutine regrid_FieldByName(src_name, src_data, dst_name, dst_data, &
                                  regridID, src_gavg, dst_gavg, &
                                  src_grid_id, dst_grid_id, &
                                  src_fld_name, dst_fld_name, &
                                  method, norm_type, line_type, renormalize_destarea, &
                                  use_src_mask, use_dst_mask, &
                                  src_mask_values, dst_mask_values, missing_value)
      character(*), intent(in)  :: src_name
      real(kind=8), intent(in)  :: src_data(:)
      character(*), intent(in)  :: dst_name
      real(kind=8), intent(out), pointer :: dst_data(:)
      integer,      intent(out), optional :: regridID
      real(ESMF_KIND_R8), intent(out), optional :: src_gavg
      real(ESMF_KIND_R8), intent(out), optional :: dst_gavg
      integer, intent(in), optional :: src_grid_id
      integer, intent(in), optional :: dst_grid_id
      character(*), intent(in), optional :: src_fld_name
      character(*), intent(in), optional :: dst_fld_name
      character(*), intent(in), optional :: method
      character(*), intent(in), optional :: norm_type
      character(*), intent(in), optional :: line_type
      logical, intent(in), optional :: renormalize_destarea
      logical, intent(in), optional :: use_src_mask
      logical, intent(in), optional :: use_dst_mask
      integer, intent(in), optional :: src_mask_values(:)
      integer, intent(in), optional :: dst_mask_values(:)
      real(kind=8), intent(in), optional :: missing_value

      !--- Local
      integer :: id_regrid
      integer :: local_src_grid_id, local_dst_grid_id, ic
      character(64) :: local_src_fld_name
      character(64) :: local_dst_fld_name
      logical :: local_renormalize_destarea
      type(cpl_vinfo_t) :: src_vinfo, dst_vinfo
      character(32) :: regrid_method
      character(32) :: local_norm_type
      character(32) :: local_line_type
      real(ESMF_KIND_R8), pointer :: tmp_data(:) => null()
      character(128) :: strng1 = " "
      character(128) :: strng2 = " "
      logical :: local_use_src_mask, local_use_dst_mask
      integer, pointer :: local_src_mask_values(:) => null()
      integer, pointer :: local_dst_mask_values(:) => null()
      real(ESMF_KIND_R8) :: local_missing_value
      real(ESMF_KIND_R8) :: local_src_gavg, local_dst_gavg

      if ( present(src_grid_id) ) then
        !--- Verify that the input src grid ID is within the valid range
        if ( src_grid_id < 1 .or. src_grid_id > n_grid_list ) then
          write(6,*)"regrid_FieldByName: Invalid src_grid_id = ",src_grid_id
          call flush(6)
          call xit("regrid_FieldByName",-1)
        endif
      endif

      if ( present(dst_grid_id) ) then
        !--- Verify that the input dst grid ID is within the valid range
        if ( dst_grid_id < 1 .or. dst_grid_id > n_grid_list ) then
          write(6,*)"regrid_FieldByName: Invalid dst_grid_id = ",dst_grid_id
          call flush(6)
          call xit("regrid_FieldByName",-2)
        endif
      endif

      !--- Determine src field name
      if ( present(src_fld_name) ) then
        !--- Allow the user to supply a name that is different from src_name
        local_src_fld_name = trim(src_fld_name)
      else
        local_src_fld_name = trim(src_name)
      endif

      !--- Determine dst field name
      if ( present(dst_fld_name) ) then
        !--- Allow the user to supply a name that is different from dst_name
        local_dst_fld_name = trim(dst_fld_name)
      else
        local_dst_fld_name = trim(dst_name)
      endif

      !--- Determine if this remap uses a mask on the source grid
      if ( present(use_src_mask) ) then
        local_use_src_mask = use_src_mask
      else
        local_use_src_mask = .false.
      endif
      if ( present(src_mask_values) ) then
        if ( size(src_mask_values) < 1 .or. size(src_mask_values) > 50 ) then
          write(6,*)"regrid_FieldByName: size(src_mask_values)=",size(src_mask_values)
          call flush(6)
          call xit("regrid_FieldByName",-3)
        endif
        allocate( local_src_mask_values(size(src_mask_values)) )
        local_src_mask_values = src_mask_values
        !--- If src_mask_values is present then always use the src mask
        local_use_src_mask = .true.
      else
        allocate( local_src_mask_values(1) )
        local_src_mask_values = 0
      endif

      !--- Determine if this remap uses a mask on the destination grid
      if ( present(use_dst_mask) ) then
        local_use_dst_mask = use_dst_mask
      else
        local_use_dst_mask = .false.
      endif
      if ( present(dst_mask_values) ) then
        if ( size(dst_mask_values) < 1 .or. size(dst_mask_values) > 50 ) then
          write(6,*)"regrid_FieldByName: size(dst_mask_values)=",size(dst_mask_values)
          call flush(6)
          call xit("regrid_FieldByName",-4)
        endif
        allocate( local_dst_mask_values(size(dst_mask_values)) )
        local_dst_mask_values = dst_mask_values
        !--- If dst_mask_value is present the always use the dst mask
        local_use_dst_mask = .true.
      else
        allocate( local_dst_mask_values(1) )
        local_dst_mask_values = 0
      endif

      !--- Determine the missing value used for unmapped destination grid cells
      if ( present(missing_value) ) then
        local_missing_value = missing_value
      else
        local_missing_value = 0.0_8
      endif

      if ( present(renormalize_destarea) ) then
        local_renormalize_destarea = renormalize_destarea
      else
        local_renormalize_destarea = .true.
      endif

      !--- Access information about the src and dst fields
      src_vinfo = find_cpl_vinfo( name=trim(adjustl(src_name)) )
      dst_vinfo = find_cpl_vinfo( name=trim(adjustl(dst_name)) )

      !--- Verify that the input src_data field is large enough
      if ( size(src_data) < src_vinfo%size ) then
        write(6,*)"regrid_FieldByName: Input source data array is too small."
        write(6,*)"  size(src_data) = ",size(src_data)," ...expecting size >= ",src_vinfo%size
        call print_cpl_vinfo_t(src_vinfo)
        call print_cpl_vinfo_t(dst_vinfo)
        call flush(6)
        call xit("regrid_FieldByName",-5)
      endif

      !--- Determine source grid ID
      if ( present(src_grid_id) ) then
        local_src_grid_id = src_grid_id
      else
        local_src_grid_id = find_grid_id( src_vinfo%grid )
        if ( local_src_grid_id < 1 .or. local_src_grid_id > n_grid_list ) then
          write(6,*)"regrid_FieldByName: Invalid src_grid_id=",local_src_grid_id," from find_grid_id"
          call flush(6)
          call xit("regrid_FieldByName",-6)
        endif
      endif

      !--- Determine destination grid ID
      if ( present(dst_grid_id) ) then
        local_dst_grid_id = dst_grid_id
      else
        local_dst_grid_id = find_grid_id( dst_vinfo%grid )
        if ( local_dst_grid_id < 1 .or. local_dst_grid_id > n_grid_list ) then
          write(6,*)"regrid_FieldByName: Invalid dst_grid_id=",local_dst_grid_id," from find_grid_id"
          call flush(6)
          call xit("regrid_FieldByName",-7)
        endif
      endif

      !--- Determine the regrid method (e.g. conservative, bilinear, ...)
      if ( present(method) ) then
        regrid_method = trim(method)
      else
        !--- Use the regrid method associated with the destination grid when
        regrid_method = trim(dst_vinfo%regrid%method)
      endif
      if ( len_trim(regrid_method) < 1 ) then
        write(6,'(3a)')'regrid_FieldByName: The regrid method for ', &
            trim(dst_name),' is missing.'
        call flush(6)
        call xit("regrid_FieldByName",-8)
      endif

      !--- Determine the norm_type (e.g. destarea, fracarea)
      if ( present(norm_type) ) then
        local_norm_type = norm_type
      else
        !--- Use the norm_type associated with the destination grid when
        local_norm_type = trim(dst_vinfo%regrid%norm_opt)
      endif

      !--- Determine the line_type (e.g. spherical, cartesian)
      if ( present(line_type) ) then
        local_line_type = line_type
      else
        !--- Default to spherical
        local_line_type = "spherical"
      endif

      !--- Never regrid a field that has the same src and dst grids
      if ( local_src_grid_id == local_dst_grid_id ) then
        !--- Verify that the sizes are also the same
        if ( src_vinfo%size /= dst_vinfo%size ) then
          write(6,'(5a,2i8)')'regrid_FieldByName: ',trim(src_name),' and ',trim(dst_name), &
            ' are on the same grid but have different sizes: ',src_vinfo%size,dst_vinfo%size
          call flush(6)
          call xit("regrid_FieldByName",-7)
        endif

        !--- Simply copy src field to dst field
        !--- Always deallocate (if allocated) then reallocate dst_data
        if ( associated(dst_data) ) deallocate(dst_data)
        allocate( dst_data(src_vinfo%size) )
        dst_data(1:src_vinfo%size) = src_data(1:src_vinfo%size)

        if ( verbose > 0 ) then
          write(6,'(6a)')'regrid_FieldByName: ',trim(src_name),' and ',trim(dst_name), &
              ' are on the same grid: ',trim(src_vinfo%grid)
          call flush(6)
        endif
        return
      endif

      !--- Find the appropriate regrid ID
      !--- Create a new regrid ID if none was found
      call find_regrid_id(id_regrid, &
                          trim(src_vinfo%grid),  &
                          trim(dst_vinfo%grid),  &
                          trim(regrid_method),   &
                          local_use_src_mask,    &
                          local_use_dst_mask,    &
                          local_src_mask_values, &
                          local_dst_mask_values, &
                          src_grid_id=local_src_grid_id, &
                          dst_grid_id=local_dst_grid_id, &
                          norm_type=local_norm_type,     &
                          line_type=local_line_type,     &
                          create=.true. )

#undef DBG_bfb
#ifdef DBG_bfb
strng1 = " "
write(strng1,'(2a)')trim(src_name),"_before_regrid_FieldById"
call dump_farray(trim(strng1),src_data(1:src_vinfo%size),idx=cpl_elapsed_time_secs)
#endif

      !--- Allocate space to hold temporary field values
      !--- This must be large enough to hold either src_data or dst_data
      allocate( tmp_data( max(src_vinfo%size,dst_vinfo%size) ) )

      tmp_data(1:src_vinfo%size) = src_data(1:src_vinfo%size)

      call regrid_FieldById(id_regrid, src_vinfo%olap, dst_vinfo%olap, tmp_data, &
                            missing_value=local_missing_value,               &
                            renormalize_destarea=local_renormalize_destarea, &
                            src_fld_name=local_src_fld_name,                 &
                            dst_fld_name=local_dst_fld_name,                 &
                            src_gavg=local_src_gavg,                         &
                            dst_gavg=local_dst_gavg)

      !--- Return the regrid ID to the caller, if requested
      if ( present(regridID) ) regridID = id_regrid

      !--- Return global averages to the caller, if requested
      if ( present(src_gavg) ) src_gavg=local_src_gavg
      if ( present(dst_gavg) ) dst_gavg=local_dst_gavg

      !--- Always deallocate (if allocated) then reallocate dst_data
      if ( associated(dst_data) ) deallocate(dst_data)
      allocate( dst_data(dst_vinfo%size) )
      dst_data(1:dst_vinfo%size) = tmp_data(1:dst_vinfo%size)

#ifdef DBG_bfb
strng1 = " "
write(strng1,'(2a)')trim(dst_name),"_after_regrid_FieldById"
call dump_farray(trim(strng1),dst_data(1:dst_vinfo%size),idx=cpl_elapsed_time_secs)
#endif

      if ( associated(tmp_data) ) nullify( tmp_data )
      if ( associated(local_src_mask_values) ) deallocate(local_src_mask_values)
      if ( associated(local_dst_mask_values) ) deallocate(local_dst_mask_values)

      if ( verbose > 0) then
        strng1 = " "
        ic = max(len_trim(src_name),len_trim(dst_name))
        ic = max(ic,10)
        strng2 = " "
        strng2 = trim(src_name)
        strng2(1:ic) = adjustr(strng2(1:ic))
        strng1 = "regrid_FieldByName: "//trim(regrid_method)//" "//trim(local_norm_type)// &
                ", SRC data "//strng2(1:ic)
        call print_var_stats(src_data, src_vinfo%size, pfx=trim(strng1), name=src_name, &
                             gavg=local_src_gavg)
        strng2 = " "
        strng2 = trim(dst_name)
        strng2(1:ic) = adjustr(strng2(1:ic))
        strng1 = "regrid_FieldByName: "//trim(regrid_method)//" "//trim(local_norm_type)// &
                ", DST data "//strng2(1:ic)
        call print_var_stats(dst_data, dst_vinfo%size, pfx=trim(strng1), name=dst_name, &
                             gavg=local_dst_gavg)
      endif

    end subroutine regrid_FieldByName

    !***************************************************************************
    !--- Regrid one field onto another field using the specifed route handle
    !***************************************************************************
    subroutine regrid_FieldById(id_regrid, src_olap, dst_olap, wrk, &
                                missing_value, renormalize_destarea, &
                                src_fld_name, dst_fld_name, &
                                src_gavg, dst_gavg)
      integer, intent(in)  :: id_regrid
      integer, intent(in)  :: src_olap
      integer, intent(in)  :: dst_olap
      real(ESMF_KIND_R8), intent(inout)  :: wrk(:)
      real(ESMF_KIND_R8), intent(in), optional :: missing_value
      logical, intent(in), optional :: renormalize_destarea
      character(*), intent(in), optional :: src_fld_name
      character(*), intent(in), optional :: dst_fld_name
      real(ESMF_KIND_R8), intent(out), optional :: src_gavg
      real(ESMF_KIND_R8), intent(out), optional :: dst_gavg

      !--- Local
      type(ESMF_Field) :: src_field, dst_field

      !--- Pointer used to access data in field objects
      real(ESMF_KIND_R8), pointer, dimension(:,:) :: a2d_ptr => null()

      !--- Temporary array
      real(ESMF_KIND_R8), pointer, dimension(:,:) :: rtmp2d => null()
      real(ESMF_KIND_R8) :: local_missing_value

      integer(ESMF_KIND_I) :: rc, lnum
      character(64) :: method = "regrid_FieldById"
      character(128) :: strng
      integer :: nx, ny, sz, idx1
      integer :: sz_in, sz_out
      integer :: src_grid_id, dst_grid_id
      logical :: local_renormalize_destarea
      real(ESMF_KIND_R8) :: local_src_gavg, local_dst_gavg, src_gsum, dst_gsum

      !--- Ensure we have a valid regrid ID
      if ( id_regrid < 1 .or. id_regrid > n_regrid_list ) then
        write(6,*)'regrid_FieldById: Invalid regrid ID = ',id_regrid
        call xit("regrid_FieldById",-1)
      endif

      if ( present(missing_value) ) then
        local_missing_value = missing_value
      else
        local_missing_value = fillvalue
      endif

      if ( present(renormalize_destarea) ) then
        local_renormalize_destarea = renormalize_destarea
      else
        local_renormalize_destarea = .true.
      endif

      !--- Assign local variable containing src and dst grid indicies
      !--- These are indicies in the global grid_list array
      src_grid_id = regrid_list(id_regrid)%src_grid_id
      dst_grid_id = regrid_list(id_regrid)%dst_grid_id

      !--- src_field will hold the incomming field values
      if ( present(src_fld_name) ) then
        src_field = ESMF_FieldCreate( grid_list(src_grid_id), ESMF_TYPEKIND_R8, rc=rc, &
                                      name=trim(src_fld_name) )
      else
        src_field = ESMF_FieldCreate( grid_list(src_grid_id), ESMF_TYPEKIND_R8, rc=rc)
      endif
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Problem creating src field", line=lnum, method=method)

      !--- dst_field will hold the interpolated field values
      if ( present(dst_fld_name) ) then
        dst_field = ESMF_FieldCreate( grid_list(dst_grid_id), ESMF_TYPEKIND_R8, rc=rc, &
                                      name=trim(dst_fld_name) )
      else
        dst_field = ESMF_FieldCreate( grid_list(dst_grid_id), ESMF_TYPEKIND_R8, rc=rc)
      endif
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Problem creating dst field", line=lnum, method=method)

      !--- Get a pointer to the 2D fortran array in src_field and
      !--- assign values to it from the incomming wrk array
      if ( associated(a2d_ptr) ) nullify( a2d_ptr )
      call ESMF_FieldGet( src_field, farrayPtr=a2d_ptr, rc=rc )
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Problem getting src_field pointer", line=lnum, method=method)

      if ( associated(a2d_ptr) ) then
        !--- Determine dimension sizes expected for this regrid operation
        nx = size(a2d_ptr,dim=1)
        ny = size(a2d_ptr,dim=2)
        sz_in = (nx + src_olap) * ny

#ifdef DBG_bfb
strng = " "
if ( present(src_fld_name) ) then
  write(strng,'(2a)')trim(src_fld_name),"_before_ESMF_FieldRegrid"
else
  write(strng,'(2a)')"srcField","_before_ESMF_FieldRegrid"
endif
call dump_farray(trim(strng),wrk(1:sz_in),idx=cpl_elapsed_time_secs)
#endif

        if ( sz_in > size(wrk) ) then
          write(6,*)"regrid_FieldById: Requesting ",sz_in, &
              " words from input data array of size ",size(wrk)
          call xit("regrid_FieldById",-2)
        endif

        if ( verbose > 2 ) then
          write(6,'(a,3i5,i8)') &
            'regrid_FieldById: Incomming field nx,ny,olap,size = ',nx,ny,src_olap,sz_in
          call flush(6)
        endif

        !--- Handle cases with extra longitudes in the incomming data
        select case (src_olap)
          case (0) !--- no overlap logitudes
            a2d_ptr = reshape( wrk(1:nx*ny), (/nx,ny/) )

          case (1) !--- 1 overlap longitude
            !--- Remove overlap longitude
            if ( associated(rtmp2d) ) deallocate(rtmp2d)
            allocate( rtmp2d(nx+src_olap,ny) )
            rtmp2d = reshape( wrk(1:(nx+src_olap)*ny), (/nx+src_olap,ny/) )
            a2d_ptr(1:nx,1:ny) = rtmp2d(1:nx,1:ny)
            deallocate( rtmp2d )

          case (2) !--- 2 overlap logitudes
            !--- Remove overlap longitudes
            if ( associated(rtmp2d) ) deallocate(rtmp2d)
            allocate( rtmp2d(nx+src_olap,ny) )
            rtmp2d = reshape( wrk(1:(nx+src_olap)*ny), (/nx+src_olap,ny/) )
            !--- resize_nemo_field will remove overlap lons
            !--- and return rtmp2d dimensioned accordingly
            call resize_nemo_field(rtmp2d)
            if ( size(rtmp2d,dim=1) == nx .and. size(rtmp2d,dim=2) == ny ) then
              a2d_ptr(1:nx,1:ny) = rtmp2d(1:nx,1:ny)
            else
              write(6,*)"regrid_FieldById: Problem removing overlap from NEMO field"
              call xit("regrid_FieldById",-3)
            endif
            deallocate( rtmp2d )

          case default  !--- Anything else is invalid
            write(6,*)"regrid_FieldById: Dimension size mismatch on source field."
            write(6,*)"    field nx, ny, olap, size = ",nx,ny,src_olap,sz_in
            call xit("regrid_FieldById",-3)

        end select

      else
        write(*,*)'Failed to get a valid src field array pointer from ESMF_FieldGet'
        call xit("regrid_FieldById",-4)
      endif

      !--- Get a pointer to the 2D fortran array in dst_field and intialize
      !--- this array with the missing value (optionally supplied by the user)
      !--- Points that do not receive a value during remapping (e.g. masked points)
      !--- will retain this value after remapping
      if ( associated(a2d_ptr) ) nullify( a2d_ptr )
      call ESMF_FieldGet( dst_field, farrayPtr=a2d_ptr, rc=rc )
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Problem getting dst_field pointer", line=lnum, method=method)

      if ( associated(a2d_ptr) ) then
        a2d_ptr = local_missing_value
      else
        write(*,*)'Failed to get a valid dst grid array pointer from ESMF_FieldGet'
        call xit("regrid_FieldById",-5)
      endif

      !--- zeroregion ---
      !---     ESMF_REGION_TOTAL  ...(default) the total regions of all DEs in dstField will
      !---                           be initialized to zero before updating the elements
      !---                           with the results of the sparse matrix multiplication.
      !---     ESMF_REGION_EMPTY  ...the elements in dstField will not be modified prior to
      !---                           the sparse matrix multiplication and results will be
      !---                           added to the incoming element values.
      !---     ESMF_REGION_SELECT ...will only zero out those elements in the destination Array
      !---                           that will be updated by the sparse matrix multiplication. 

      !--- termorderflag ---
      !--- Specifies the order of the source side terms in all of the destination sums.
      !--- The termorderflag only affects the order of terms during the execution of the
      !--- RouteHandle. The default is ESMF_TERMORDER_FREE, allowing maximum flexibility
      !--- in the order of terms for optimum performance.
      !---     ESMF_TERMORDER_SRCSEQ ...The source terms are in strict ascending order
      !---                              according to their source sequence index.
      !---     ESMF_TERMORDER_SRCPET ...The source terms are first ordered according to
      !---                              their distribution across the source side PETs:
      !---                              for each destination PET the source PET order starts
      !---                              with the localPet and decrements from there, modulo
      !---                              petCount, until all petCount PETs are accounted for.
      !---                              The term order within each source PET is given by the
      !---                              source term sequence index.
      !---     ESMF_TERMORDER_FREE   ...There is no prescribed term order. The source terms
      !---                              may be summed in any order that optimizes performance.

      !--- checkflag ---
      !--- If checkflag is set to .TRUE. the input Array pair will be checked for
      !--- consistency with the precomputed operation provided by routehandle.
      !--- If checkflag is set to .FALSE. (default) only a very basic input check
      !--- will be performed, leaving many inconsistencies undetected.
      !--- Set checkflag to .FALSE. to achieve highest performance.

      !--- Regrid this field
      call ESMF_FieldRegrid( routehandle = regrid_rh_list(id_regrid), &
                                srcField = src_field, &
                                dstField = dst_field, &
                              zeroregion = ESMF_REGION_SELECT,    &
                           termorderflag = ESMF_TERMORDER_SRCSEQ, &
                               checkflag = T_esmf,                &
                                      rc = rc )

      !--- a2d_ptr should still be pointing at the fortran array in dst_field
      if ( associated(a2d_ptr) ) then
        nx = size(a2d_ptr,dim=1)
        ny = size(a2d_ptr,dim=2)
        sz_out = (nx + dst_olap) * ny

#ifdef DBG_bfb
strng = " "
if ( present(dst_fld_name) ) then
  write(strng,'(2a)')trim(dst_fld_name),"_after_ESMF_FieldRegrid"
else
  write(strng,'(2a)')"dstField","_after_ESMF_FieldRegrid"
endif
if (associated(rtmp2d)) deallocate(rtmp2d)
allocate( rtmp2d(nx,ny) )
rtmp2d(1:nx,1:ny) = a2d_ptr(1:nx,1:ny)
call dump_farray(trim(strng),rtmp2d,idx=cpl_elapsed_time_secs)
deallocate(rtmp2d)
#endif

        !--- For weights generated using destination area normalization, the values of the
        !--- destination cells which extend partway outside the unmasked source field
        !--- are decreased by the fraction they extend outside.
        !--- To correct these values, the destination field resulting from the the call to
        !--- ESMF_FieldRegrid must be divided by the destination fraction (ie dst_frac which
        !--- is returned from the call to ESMF_FieldRegridStore)
        if ( trim(adjustl(regrid_list(id_regrid)%regrid_method)) .eq. "conservative" .and. &
             trim(adjustl(regrid_list(id_regrid)%norm_type)) .eq. "destarea" ) then
          if ( local_renormalize_destarea ) then
            !--- Renormalize destination values by destination fraction
            if ( size(regrid_list(id_regrid)%dst_frac,dim=1) /= nx .or. &
                 size(regrid_list(id_regrid)%dst_frac,dim=2) /= ny ) then
              write(6,*)"regrid_FieldById: Size mismatch between current field and dst_frac"
              write(6,*)"    shape(field   ) = ",shape(a2d_ptr)
              write(6,*)"    shape(dst_frac) = ",shape(regrid_list(id_regrid)%dst_frac)
              call xit("regrid_FieldById",-6)
            endif
            write(6,*)"regrid_FieldById: Regrid ID=",id_regrid, &
                      "  Renormalize with destination fraction."
            call flush(6)
            !--- For each grid cell: if dst fraction is non-zero then divide by dst fraction
            where ( regrid_list(id_regrid)%dst_frac .ne. 0.0 ) &
              a2d_ptr = a2d_ptr/regrid_list(id_regrid)%dst_frac
          else
            write(6,*)"regrid_FieldById: Regrid ID=",id_regrid, &
                      "  Do not renormalize with destination fraction."
            call flush(6)
          endif
        endif

        if ( sz_out > size(wrk) ) then
          write(6,*)"regrid_FieldById: Writing ",sz_out, &
              " words to output data array of size ",size(wrk)
          call xit("regrid_FieldById",-6)
        endif

        if ( verbose > 2 ) then
          write(6,'(a,3i5,i8)') &
            'regrid_FieldById:  Outgoing field nx,ny,olap,size = ',nx,ny,dst_olap,sz_out
          call flush(6)
        endif

        !--- Zero out the wrk array before reassigning
        wrk(:) = 0.0

        !--- Handle cases with extra longitudes in the outgoing data
        select case (dst_olap)
          case (0) !--- no overlap logitudes
            wrk(1:sz_out) = reshape( a2d_ptr, (/sz_out/) )

          case (1) !--- 1 overlap logitude
            !--- Add overlap longitudes to the outgoing data
            if ( associated(rtmp2d) ) deallocate(rtmp2d)
            allocate( rtmp2d(nx+dst_olap, ny) )
            rtmp2d(1:nx,1:ny) = a2d_ptr(1:nx,1:ny)
            rtmp2d(nx+1,1:ny) = rtmp2d(1,1:ny)
            wrk(1:sz_out) = reshape( rtmp2d, (/sz_out/) )
            deallocate( rtmp2d )

          case (2) !--- 2 overlap logitudes
            !--- Add overlap longitudes
            if ( associated(rtmp2d) ) deallocate(rtmp2d)
            allocate( rtmp2d(nx+dst_olap,ny) )
            !--- Assign the interior points (without extra longitudes)
            rtmp2d(2:nx+1,1:ny) = a2d_ptr(1:nx,1:ny)
            !--- Then assign the extra longitudes
            rtmp2d(1,1:ny)    = a2d_ptr(nx,1:ny)
            rtmp2d(nx+2,1:ny) = a2d_ptr(1,1:ny)
            wrk(1:sz_out) = reshape( rtmp2d, (/sz_out/) )
            deallocate( rtmp2d )

          case default  !--- Anything else is invalid
            write(6,*)"regrid_FieldById: Dimension size mismatch on destination field."
            write(6,*)"    field nx, ny, olap, size = ",nx,ny,dst_olap,sz_out
            call xit("regrid_FieldById",-7)

        end select

      else
        write(*,*)'Failed to get a valid dst grid array pointer from ESMF_FieldGet'
        call xit("regrid_FieldById",-8)
      endif

      !--- Calculate and print global averages for source and destinstion fields
      if ( present(src_gavg) ) then
        src_gavg = calc_global_integral(src_field, id_regrid, "src", 2, src_gsum)
      else
        local_src_gavg = calc_global_integral(src_field, id_regrid, "src", 2, src_gsum)
      endif
      if ( present(dst_gavg) ) then
        dst_gavg = calc_global_integral(dst_field, id_regrid, "dst", 2, dst_gsum)
      else
        local_dst_gavg = calc_global_integral(dst_field, id_regrid, "dst", 2, dst_gsum)
      endif

      !--- Release resources used by src_field and dst_field
      if ( associated(a2d_ptr) ) nullify( a2d_ptr )

      call ESMF_FieldDestroy(src_field, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to destroy src_field", line=lnum, method=method)

      call ESMF_FieldDestroy(dst_field, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to destroy dst_field", line=lnum, method=method)

    end subroutine regrid_FieldById

    !***************************************************************************
    !--- Calculate the global integral and average for field given the regrid ID
    !--- Note that ESMF_Field contains a pointer to an ESMF_FieldType class so
    !--- that passing an ESMF_Field effectively passes a pointer (ie not data)
    !***************************************************************************
    function calc_global_integral(field_in, regridID, src_or_dst, mode, gsum) result(gavg)
      !--- A field containing the data for which the global average is required
      type(ESMF_Field), intent(in) :: field_in

      !--- regridID is used to find associated area and frac fields
      integer, intent(in)          :: regridID

      !--- This string determines if the source or the destination area/frac
      !--- associated with this regridID is to be used
      character(*), intent(in)     :: src_or_dst

      !--- mode determines whether or not this routine will calculate the avg
      !---   mode < 1  ...do nothing
      !---   mode == 1 ...calculate and return gsum and gavg
      !---   mode > 1  ...calculate and return gsum and gavg and print gsum and gavg
      integer, intent(in)          :: mode

      !--- Return the global integral as gsum
      real(ESMF_KIND_R8) :: gsum

      !--- Return the global average as gavg
      real(ESMF_KIND_R8) :: gavg

      !--- Local
      real(ESMF_KIND_R8), pointer :: data_ptr(:,:) => null()
      real(ESMF_KIND_R8), pointer :: frac_ptr(:,:) => null()
      real(ESMF_KIND_R8), pointer :: area_ptr(:,:) => null()
      integer(ESMF_KIND_I) :: rc, lnum
      character(64) :: method = "calc_global_integral"
      logical :: use_src, use_dst
      integer :: nx, ny, ic
      character(32) :: fld_name
      character(64) :: strng
      real(kind=8) :: garea

      !--- Return zero when avg is not calculated
      gavg = 0.0_8
      gsum = 0.0_8

      !--- Do nothing unless mode > 0
      if ( mode < 1 ) return

      if ( regridID < 1 .or. regridID > n_regrid_list ) then
        write(6,*)"calc_global_integral: regridID=",regridID," is out of range."
        call xit("calc_global_integral",-1)
      endif

      use_src = .false.
      use_dst = .false.
      if ( trim(adjustl(src_or_dst)) .eq. "src" ) then
        use_src = .true.
      else if ( trim(adjustl(src_or_dst)) .eq. "dst" ) then
        use_dst = .true.
      else
        write(6,*)"calc_global_integral: src_or_dst must be either src or dst."
        call xit("calc_global_integral",-2)
      endif

      !--- Get a pointer to the field data
      call ESMF_FieldGet( field_in, farrayPtr=data_ptr, rc=rc )
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Problem getting field data pointer", line=lnum, method=method)
      if ( .not. associated(data_ptr) ) then
        write(6,*)"calc_global_integral: Unable to get pointer to field data."
        call xit("calc_global_integral",-3)
      endif

      !--- Determine dimension sizes from the data pointer
      nx = size(data_ptr, dim=1)
      ny = size(data_ptr, dim=2)

      !--- Get the name of this field
      fld_name = " "
      call ESMF_FieldGet( field_in, name=fld_name, rc=rc )
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Problem getting field name", line=lnum, method=method)

      !--- Get pointers to the area and frac data
      if ( use_src ) then
        if ( .not. associated(regrid_list(regridID)%src_frac) ) then
          write(6,*)"calc_global_integral: src_frac is not associated for regridID=",regridID
          call xit("calc_global_integral",-4)
        endif
        if ( .not. associated(regrid_list(regridID)%src_area) ) then
          write(6,*)"calc_global_integral: src_area is not associated for regridID=",regridID
          call xit("calc_global_integral",-5)
        endif
        frac_ptr => regrid_list(regridID)%src_frac
        area_ptr => regrid_list(regridID)%src_area
      endif

      if ( use_dst ) then
        if ( .not. associated(regrid_list(regridID)%dst_frac) ) then
          write(6,*)"calc_global_integral: dst_frac is not associated for regridID=",regridID
          call xit("calc_global_integral",-6)
        endif
        if ( .not. associated(regrid_list(regridID)%dst_area) ) then
          write(6,*)"calc_global_integral: dst_area is not associated for regridID=",regridID
          call xit("calc_global_integral",-7)
        endif
        frac_ptr => regrid_list(regridID)%dst_frac
        area_ptr => regrid_list(regridID)%dst_area
      endif

      if ( size(frac_ptr, dim=1) == nx .and. size(frac_ptr, dim=2) == ny .and. &
           size(area_ptr, dim=1) == nx .and. size(area_ptr, dim=2) == ny ) then
        !--- Calculate the global average
        gsum  = sum(data_ptr * frac_ptr * area_ptr)
        garea = sum(frac_ptr * area_ptr)
        if ( garea == 0.0_8 ) then
          gavg  = 0.0_8
        else
          gavg  = gsum/garea
        endif
      else
        write(6,*)"calc_global_integral: Shape mismatch."
        write(6,*)" shape(data) = ",nx,ny
        write(6,*)" shape(frac) = ",size(frac_ptr,dim=1),size(frac_ptr,dim=2)
        write(6,*)" shape(area) = ",size(area_ptr,dim=1),size(area_ptr,dim=2)
        call xit("calc_global_integral",-8)
      endif

      nullify(data_ptr)
      nullify(frac_ptr)
      nullify(area_ptr)

      if ( mode > 1 ) then
        !--- Print the average here
        strng = " "
        strng = trim(fld_name)
        strng = adjustr(strng)
        !--- The field name has a maximum length of 32 characters
        ic = len(strng) - 31
        write(6,'(a,a32,3(a,g23.15))') &
            "global integral: ",strng(ic:len(strng))," = ",gsum,"  avg = ",gavg,"  remap area = ",garea
        call flush(6)
      endif

    end function calc_global_integral

    !***************************************************************************
    !--- Find the grid ID given a grid name and, optionally, a mask
    !--- If no initialized grid is found then return grid ID = 0 unless the
    !--- grid name is one of atm, ocn, canom (in which case we abort)
    !***************************************************************************
    function find_grid_id(grid_name, imask) result(grid_id)
      !--- The name associated with this grid
      character(*) :: grid_name

      !--- The integer mask associated with this grid
      integer, intent(in), optional :: imask(:,:)

      !--- Return the grid ID
      integer :: grid_id

      !--- Local
      integer :: idx
      character(128) :: file_name, sfx

      grid_id = 0
      select case ( trim(adjustl(grid_name)) )
        case ("atm")
          if ( id_atm_grid < 0 ) then
            write(6,'(a)')"find_grid_id: Uninitialized atm grid."
            call flush(6)
            call xit("find_grid_id",-1)
          endif
          grid_id = id_atm_grid

        case ("ocn")
          if ( id_ocn_grid < 0 ) then
            write(6,'(a)')"find_grid_id: Uninitialized ocn grid."
            call flush(6)
            call xit("find_grid_id",-2)
          endif
          grid_id = id_ocn_grid

        case ("canom")
          if ( id_canom4_grid < 0 ) then
            write(6,'(a)')"find_grid_id: Uninitialized canom grid."
            call flush(6)
            call xit("find_grid_id",-3)
          endif
          grid_id = id_canom4_grid

        case default
          !--- Is this a file name
          !--- A missing grid does not result in an abort in this case
          if ( associated(cpl_grid) ) then
            do idx=1,size(cpl_grid)
              !--- Ignore grids without an associated file name
              if ( len_trim(cpl_grid(idx)%file_name) < 1 ) cycle
              file_name = " "
              sfx = " "
              file_name = trim( grid_desc_file_name(idx, sfx=sfx) )
              if ( trim(adjustl(file_name)) .eq. trim(adjustl(grid_name)) ) then
                !--- This is a file name associated with an existing grid
                if ( present(imask) ) then
                  if ( .not. associated(cpl_grid(idx)%mask) ) then
                    !--- If the internal mask is not available then accept this match
                    write(6,*)"find_grid_id: Request to match imask but imask is not" &
                              //" defined for grid ",trim(file_name)
                    write(6,*)"    Ignoring match request and using grid ", &
                              trim(file_name)
                    grid_id = idx
                    exit
                  endif
                  !--- The input mask and the mask stored in cpl_grid must match
                  if ( size(imask,dim=1) == size(cpl_grid(idx)%mask,dim=1) .and. &
                       size(imask,dim=2) == size(cpl_grid(idx)%mask,dim=2) ) then
                    if ( all(imask == cpl_grid(idx)%mask) ) then
                      !--- imask is a match
                      grid_id = idx
                      exit
                    endif
                  else
                    !--- The shape of imask does not match the
                    !--- shape of the mask in cpl_grid
                    write(6,*)"find_grid_id: shape mismatch for imask on grid ",trim(file_name)
                    write(6,*)"    found shape(imask) = ",shape(imask)
                    write(6,*)" expected shape(imask) = ",shape(cpl_grid(idx)%mask)
                  endif
                else
                  !--- Accept this match
                  grid_id = idx
                  exit
                endif
              endif
            enddo
          endif
      end select

    end function find_grid_id

    !***************************************************************************
    !--- Given a grid ID return the associated grid description file name
    !***************************************************************************
    function grid_desc_file_name(grid_id, sfx, fid) result(file_name)
      integer, intent(in) :: grid_id
      character(*), intent(out), optional :: sfx
      integer,      intent(out), optional :: fid

      !--- The return value is the file name without any suffix
      character(128) :: file_name

      !--- Local
      integer :: c1, c2, idx

      if ( grid_id < 1 .or. grid_id > n_grid_list ) then
        write(6,*)"grid_desc_file_name: grid_id = ",grid_id," is out of range."
        call xit("grid_desc_file_name",-1)
      endif

      !--- Get the file name from cpl_grid
      file_name = " "
      file_name = trim(cpl_grid(grid_id)%file_name)
      c2 = len_trim(file_name)
      if ( c2 < 1 ) then
        !--- If file_name is blank then simply return blank
        if ( present(sfx) ) sfx = " "
        if ( present(fid) ) fid = 0
        return
      endif

      if ( file_name(c2:c2) .eq. ")" ) then
        !--- Assume there is a suffix of the form (...) on this file name
        !--- which is not part of the file name found on disk
        c1 = index(file_name,"(",back=.true.)
        if ( c1 > 0 ) then
          !--- Strip the suffix from file_name
          if ( present(sfx) ) sfx = file_name(c1:c2)
          if ( present(fid) ) then
            if ( c2-c1 < 2 ) then
              !--- There are no characters between the parentheses in sfx
              write(6,*)"grid_desc_file_name: Empty suffix on ",trim(file_name)
              call flush(6)
              call xit("grid_desc_file_name",-2)
            endif
            if (verify(file_name(c1+1:c2-1),'0123456789') .ne. 0) then
              !--- Verify that all characters between the parentheses are 0-9
              write(6,*)"grid_desc_file_name: Invalid suffix on ",trim(file_name)
              call flush(6)
              call xit("grid_desc_file_name",-3)
            endif
            !--- Assign the file ID
            read(file_name(c1+1:c2-1),*) fid
          endif
          file_name = file_name(1:c1-1)
        else
          !--- Unmatched parentheses
          write(6,*)"grid_desc_file_name: Unmatched parentheses on ",trim(file_name)
          call flush(6)
          call xit("grid_desc_file_name",-4)
        endif
      else
        !--- This file name does not contain a suffix of the form (...)
        if ( present(sfx) ) sfx = " "
        if ( present(fid) ) fid = 0
      endif

    end function grid_desc_file_name

    !***************************************************************************
    !--- Create grid objects from SCRIP format grid description files
    !***************************************************************************
    subroutine init_grids()

      !--- Local
      integer(ESMF_KIND_I) :: rc
      integer(ESMF_KIND_I) :: positionVector(2)
      integer(ESMF_KIND_I) :: nx, ny
      type(ESMF_DistGridConnection), pointer :: connectionList(:) => null()
      type(ESMF_DistGrid) :: distGrid
      logical :: exists
      character(128) :: grid_title

      !--- Allocate the grid list
      !--- The first 4 elements in this list are reserved
      n_grid_list = 4
      if ( .not. associated(grid_list) ) then
        !--- No grids have been defined yet
        allocate( grid_list(max_grid_list) )
        if ( associated(cpl_grid) ) nullify(cpl_grid)
        allocate( cpl_grid(max_grid_list) )
      endif

      !--- Read the default atmosphere grid description file if present
      id_atm_grid = -1
      if ( len_trim(atm_grid_desc) > 0 ) then
        id_atm_grid = 1
        call define_grid_from_file(id_atm_grid, trim(atm_grid_desc))

        if (verbose > 0) then
          write(6,'(2a)')'init_grids: Initialized atm grid from file ',trim(atm_grid_desc)
          call flush(6)
        endif

        if (verbose > 1) then
          call print_grid_info( grid_list(id_atm_grid) )
        endif
      else
        !--- Define the atm grid using coordinate info passed from the AGCM at runtime
        id_atm_grid = 1
        call define_grid_from_agcm_coords(id_atm_grid)

        if (verbose > 0) then
          write(6,'(2a)')'init_grids: Initialized atm grid from coord info passed from the AGCM'
          call flush(6)
        endif
      endif

      !--- Read the default ocean grid description file if present
      id_ocn_grid = -1
      if ( len_trim(ocn_grid_desc) > 0 ) then
        id_ocn_grid = 2
        call define_grid_from_file(id_ocn_grid, trim(ocn_grid_desc))

        if (verbose > 0) then
          write(6,'(2a)')'init_grids: Initialized ocn grid from file ',trim(ocn_grid_desc)
          call flush(6)
        endif

        if (verbose > 1) then
          call print_grid_info( grid_list(id_ocn_grid) )
        endif
      else
        !--- Define the ocean grid using coordinate info passed from NEMO at runtime
        id_ocn_grid = 2
        call define_grid_from_nemo_coords(id_ocn_grid)

        if (verbose > 0) then
          write(6,'(2a)')'init_grids: Initialized ocean grid from coord info passed from NEMO'
          call flush(6)
        endif
      endif

      if (verbose > 1) then
        call print_grid_info( grid_list(id_atm_grid) )
        call print_grid_info( grid_list(id_ocn_grid) )
      endif

      !--- Read the default ice grid desciption file if present
      id_ice_grid = -1
      if ( len_trim(ice_grid_desc) > 0 ) then
        id_ice_grid = 3
        call define_grid_from_file(id_ice_grid, trim(ice_grid_desc))

        if (verbose > 0) then
          write(6,'(2a)')'init_grids: Initialized ice grid from file ',trim(ice_grid_desc)
          call flush(6)
        endif

        if (verbose > 1) then
          call print_grid_info( grid_list(id_ice_grid) )
        endif
      endif

      !--- Read the default canom grid desciption file if present
      id_canom4_grid = -1
      if ( len_trim(canom4_grid_desc) > 0 ) then
        id_canom4_grid = 4
        call define_grid_from_file(id_canom4_grid, trim(canom4_grid_desc))

        if (verbose > 0) then
          write(6,'(2a)')'init_grids: Initialized canom grid from file ',trim(canom4_grid_desc)
          call flush(6)
        endif

        if (verbose > 1) then
          call print_grid_info( grid_list(id_canom4_grid) )
        endif
      endif

    end subroutine init_grids

    !***************************************************************************
    !--- Add a grid object using info from a grid description file
    !***************************************************************************
    subroutine add_grid(grid_id, grid_desc_file, imask, overlap, cartesian, &
                        grid_dims, grid_edges, grid_size)
      !--- The index in grid_list and cpl_grid of this grid
      integer, intent(out) :: grid_id

      !--- The name of the SCRIP format grid description file
      character(*), intent(in) :: grid_desc_file

      !--- The integer mask associated with this grid
      !--- This will overwrite the mask read from the grid description file
      integer, intent(in), optional :: imask(:,:)

      !--- Include overlap longitude (if found in file) in the new grid object
      logical, intent(in), optional :: overlap

      !--- Create a cartesian grid with no periodic dimensions
      logical, intent(in), optional :: cartesian

      integer, intent(out), optional :: grid_dims(2)
      integer, intent(out), optional :: grid_edges(2,2)
      integer, intent(out), optional :: grid_size

      !--- Local
      character(64) :: method = "add_grid"
      integer(ESMF_KIND_I) :: lnum
      integer(ESMF_KIND_I) :: rc
      integer :: idx
      type(ESMF_DistGrid) :: distgrid
      character(128) :: grid_name
      integer :: local_grid_dims(2)
      integer :: local_grid_edges(2,2)
      integer :: local_grid_size
      logical :: cart

      if ( present(cartesian) ) then
        cart = cartesian
      else
        cart = .false.
      endif

      if ( .not. associated(grid_list) ) then
        !--- No grids have been defined yet
        n_grid_list = 1
        grid_id = 1
        allocate( grid_list(max_grid_list) )
        if ( associated(cpl_grid) ) nullify(cpl_grid)
        allocate( cpl_grid(max_grid_list) )

      else
        !--- Determine if this grid is already defined
        if ( present(imask) ) then
          grid_id = find_grid_id(grid_desc_file, imask=imask)
        else
          grid_id = find_grid_id(grid_desc_file)
        endif

        !--- If it is defined then return the grid ID along with any optional output
        if ( grid_id > 0 ) then
          if ( present(grid_dims) ) grid_dims(1:2) = &
                                    (/ cpl_grid(grid_id)%nx, cpl_grid(grid_id)%ny /)
          if ( present(grid_edges) ) grid_edges = cpl_grid(grid_id)%edges
          if ( present(grid_size) ) grid_size = cpl_grid(grid_id)%size
          return
        endif

        !--- Increment the total number of grid definitions in grid_list
        n_grid_list = n_grid_list + 1

        !--- Set grid_id to the next unassigned element of grid_list
        grid_id = n_grid_list

      endif

      if ( grid_id < 1 .or. grid_id > max_grid_list ) then
        write(6,*)"add_grid: grid ID = ",grid_id," is out of range."
        call xit("add_grid",-1)
      endif

      if ( present(imask) .and. present(overlap) ) then
        call define_grid_from_file(grid_id, grid_desc_file, &
                                   imask=imask, &
                                   overlap=overlap, &
                                   grid_dims=local_grid_dims, &
                                   grid_edges=local_grid_edges, &
                                   grid_size=local_grid_size, &
                                   cartesian=cart)
      else if ( present(imask) ) then
        call define_grid_from_file(grid_id, grid_desc_file, &
                                   imask=imask, &
                                   grid_dims=local_grid_dims, &
                                   grid_edges=local_grid_edges, &
                                   grid_size=local_grid_size, &
                                   cartesian=cart)
      else if ( present(overlap) ) then
        call define_grid_from_file(grid_id, grid_desc_file, &
                                   overlap=overlap, &
                                   grid_dims=local_grid_dims, &
                                   grid_edges=local_grid_edges, &
                                   grid_size=local_grid_size, &
                                   cartesian=cart)
      else
        call define_grid_from_file(grid_id, grid_desc_file, &
                                   grid_dims=local_grid_dims, &
                                   grid_edges=local_grid_edges, &
                                   grid_size=local_grid_size, &
                                   cartesian=cart)
      endif

      if ( present(grid_dims) )  grid_dims  = local_grid_dims
      if ( present(grid_edges) ) grid_edges = local_grid_edges
      if ( present(grid_size) )  grid_size  = local_grid_size

      write(6,*)"add_grid: New grid ID = ",grid_id, &
                "   name = ",trim(cpl_grid(grid_id)%file_name)
      if (verbose > 1) then
        call print_grid_info( grid_list(grid_id) )
      endif

    end subroutine add_grid

    !***************************************************************************
    !--- Define a grid object using data read from a grid description file
    !***************************************************************************
    subroutine define_grid_from_file(grid_id, file_name, imask, use_mask, overlap, &
                                     grid_dims, grid_edges, grid_size, cartesian)
      !--- The index in grid_list and cpl_grid of this grid
      integer, intent(in) :: grid_id

      !--- The name of the SCRIP format grid description file
      character(*), intent(in) :: file_name

      !--- The integer mask associated with this grid
      !--- This will overwrite the mask read from the grid description file
      integer, intent(in), optional :: imask(:,:)

      !--- use_mask = False means imask is set to 1 everywhere
      !--- use_mask = True  means imask is what was read from the file
      logical, intent(in), optional :: use_mask

      !--- Include overlap longitude (if found in file) in the new grid object
      logical, intent(in), optional :: overlap

      integer, intent(out), optional :: grid_dims(2)
      integer, intent(out), optional :: grid_edges(2,2)
      integer, intent(out), optional :: grid_size

      !--- Create a cartesian grid with no periodic dimensions
      logical, intent(in), optional :: cartesian

      !--- Local
      character(64) :: method = "define_grid_from_file"
      integer(ESMF_KIND_I) :: lnum
      integer(ESMF_KIND_I) :: rc
      integer(ESMF_KIND_I) :: regDecomp(2)
      integer(ESMF_KIND_I) :: positionVector(2)
      integer(ESMF_KIND_I) :: nx, ny
      integer(ESMF_KIND_I) :: nxc, nyc
      integer(ESMF_KIND_I) :: gridEdgeLWidth(2), gridEdgeUWidth(2)
      type(ESMF_DistGridConnection) :: connectionList(1)
      type(ESMF_DistGrid) :: distGrid
      logical :: exists, cart, add_area_grid_item, file_was_accessed=.false.
      integer :: idx1, idx2, olapx, fid, fidmax
      character(128) :: sfx, strng

      !--- Info to be read from the grid description file
      character(128)         :: grid_title
      integer                :: local_grid_dims(2)
      integer                :: local_grid_edges(2,2)
      integer                :: local_grid_size
      integer,       pointer :: local_imask(:,:) => null()
      real(kind=8),  pointer :: mid_lon(:,:) => null()
      real(kind=8),  pointer :: mid_lat(:,:) => null()
      real(kind=8),  pointer :: corn_lon(:,:) => null()
      real(kind=8),  pointer :: corn_lat(:,:) => null()
      real(kind=8),  pointer :: mid_lon_1d(:) => null()
      real(kind=8),  pointer :: mid_lat_1d(:) => null()
      real(kind=8),  pointer :: corn_lon_1d(:) => null()
      real(kind=8),  pointer :: corn_lat_1d(:) => null()
      integer :: grid_file_type

      !--- Internal info associated with this grid
      real(ESMF_KIND_R8), pointer    :: grid_mid_lon(:,:) => null()
      real(ESMF_KIND_R8), pointer    :: grid_mid_lat(:,:) => null()
      real(ESMF_KIND_R8), pointer    :: grid_corn_lon(:,:) => null()
      real(ESMF_KIND_R8), pointer    :: grid_corn_lat(:,:) => null()
      integer(ESMF_KIND_I4), pointer :: grid_mask(:,:) => null()

      real(ESMF_KIND_R8), pointer    :: grid_mid_lon_1d(:) => null()
      real(ESMF_KIND_R8), pointer    :: grid_mid_lat_1d(:) => null()
      real(ESMF_KIND_R8), pointer    :: grid_corn_lon_1d(:) => null()
      real(ESMF_KIND_R8), pointer    :: grid_corn_lat_1d(:) => null()

      !--- Grid area related
      type(ESMF_Field) :: area_field
      real(ESMF_KIND_R8), pointer :: area_data1(:,:) => null()
      real(ESMF_KIND_R8), pointer :: area_data2(:,:) => null()
      real(ESMF_KIND_R8), pointer :: area_data3(:,:) => null()
      real(ESMF_KIND_R8), pointer :: mask_data(:,:) => null()
      real(ESMF_KIND_R8) :: garea, rad2

      !--- Ensure that this file exists
      if ( .not. locate_file(trim(file_name), remote=file_was_accessed) ) then
        write(6,*)'define_grid_from_file: Missing grid description file ',trim(file_name)
        call xit("define_grid_from_file",-1)
      endif

      if ( present(cartesian) ) then
        cart = cartesian
      else
        cart = .false.
      endif

      !--- add_area_grid_item flags the addition of grid cell area to the current grid
      add_area_grid_item = .true.

      !--- Define regular decomposition counts for each dimension
      regDecomp(1) = cpl_PetCount
      regDecomp(2) = 1

      !--- Identify the type of grid_description file that is to be read
      grid_file_type = grid_desc_file_type(file_name)

      !--- Read the grid description file
      if ( present(use_mask) .and. present(overlap) ) then
        call read_scrip_grid_file(file_name, title=grid_title, grid_size=local_grid_size, &
                                  grid_dims=local_grid_dims, imask=local_imask, &
                                  mid_lon=mid_lon, mid_lat=mid_lat, &
                                  corn_lon=corn_lon, corn_lat=corn_lat, &
                                  grid_edges=local_grid_edges, &
                                  use_mask=use_mask, overlap=overlap)
      else if ( present(use_mask) ) then
        call read_scrip_grid_file(file_name, title=grid_title, grid_size=local_grid_size, &
                                  grid_dims=local_grid_dims, imask=local_imask, &
                                  mid_lon=mid_lon, mid_lat=mid_lat, &
                                  corn_lon=corn_lon, corn_lat=corn_lat, &
                                  grid_edges=local_grid_edges, &
                                  use_mask=use_mask)
      else if ( present(overlap) ) then
        call read_scrip_grid_file(file_name, title=grid_title, grid_size=local_grid_size, &
                                  grid_dims=local_grid_dims, imask=local_imask, &
                                  mid_lon=mid_lon, mid_lat=mid_lat, &
                                  corn_lon=corn_lon, corn_lat=corn_lat, &
                                  grid_edges=local_grid_edges, &
                                  overlap=overlap)
      else
        select case (grid_file_type)
          case (1) !--- A SCRIP format file
            call read_scrip_grid_file(file_name, title=grid_title, grid_size=local_grid_size, &
                                  grid_dims=local_grid_dims, imask=local_imask, &
                                  mid_lon=mid_lon, mid_lat=mid_lat, &
                                  mid_lon_1d=mid_lon_1d, mid_lat_1d=mid_lat_1d, &
                                  corn_lon=corn_lon, corn_lat=corn_lat, &
                                  corn_lon_1d=corn_lon_1d, corn_lat_1d=corn_lat_1d, &
                                  grid_edges=local_grid_edges)
          case (2) !--- A NEMO "mesh_mask" format file
            call read_NEMO_mesh_mask_file(file_name, title=grid_title, grid_size=local_grid_size, &
                                  grid_dims=local_grid_dims, imask=local_imask, &
                                  mid_lon=mid_lon, mid_lat=mid_lat, &
                                  corn_lon=corn_lon, corn_lat=corn_lat, &
                                  grid_edges=local_grid_edges)
          case default
            write(6,*)'define_grid_from_file: Invalid grid_file_type = ',grid_file_type
            call xit("define_grid_from_file",-1)
        end select
      endif

      if ( present(grid_dims) )  grid_dims  = local_grid_dims
      if ( present(grid_edges) ) grid_edges = local_grid_edges
      if ( present(grid_size) )  grid_size  = local_grid_size
      nx = local_grid_dims(1)
      ny = local_grid_dims(2)
      olapx = local_grid_edges(2,1)
      write(6,'(a,i4,3i8)')"define_grid_from_file: olapx,nx,ny,size = ", &
                            olapx,nx,ny,local_grid_size
      call flush(6)

      if ( present(imask) .and. associated(local_imask) ) then
        !--- Overwrite the mask read from the file with a mask supplied
        !--- via the parameter list for this routine

        !--- Ensure these masks have the same shape
        if ( size(imask,dim=1) /= size(local_imask,dim=1) .or. &
             size(imask,dim=2) /= size(local_imask,dim=2) ) then
          write(6,*)"define_grid_from_file: Shape mismatch for mask areas"
          write(6,*)"   found nx,ny = ",size(imask,dim=1),size(imask,dim=2)
          write(6,*)"expected nx,ny = ",size(local_imask,dim=1),size(local_imask,dim=2)
          call xit("define_grid_from_file",-1)
        endif

        !--- Overwrite with the new mask values
        local_imask = imask
      endif

      if ( verbose > 10 ) then
        !--- Print some coordinate info
        write(6,'(a,2i8)')"define_grid_from_file:    shape(imask) = ",shape(local_imask)
        write(6,'(a,2i8)')"define_grid_from_file:  shape(mid_lon) = ",shape(mid_lon)
        write(6,'(a,2i8)')"define_grid_from_file:  shape(mid_lat) = ",shape(mid_lat)
        write(6,'(a,2i8)')"define_grid_from_file: shape(corn_lon) = ",shape(corn_lon)
        write(6,'(a,2i8)')"define_grid_from_file: shape(corn_lat) = ",shape(corn_lat)
        do idx2=1,ny
          if ( idx2 > 2 .and. idx2 < ny-1 ) cycle
          do idx1=1,nx
            if ( idx1 > 2 .and. idx1 < nx-1 ) cycle
            if ( idx1 == nx ) then
              write(6,'(2i6,5(1x,a,2g16.8,a))')idx1,idx2, &
                      ": (",mid_lon(idx1,idx2),mid_lat(idx1,idx2),")",   &
                      "corners 1 and 4 (",corn_lon(idx1,idx2),corn_lat(idx1,idx2),")", &
                      "(",corn_lon(idx1,idx2+1),corn_lat(idx1,idx2+1),")"
            else
              write(6,'(2i6,5(1x,a,2g16.8,a))')idx1,idx2, &
                      ": (",mid_lon(idx1,idx2),mid_lat(idx1,idx2),")",   &
                      "corners (",corn_lon(idx1,idx2),corn_lat(idx1,idx2),")", &
                      "(",corn_lon(idx1+1,idx2),corn_lat(idx1+1,idx2),")",   &
                      "(",corn_lon(idx1+1,idx2+1),corn_lat(idx1+1,idx2+1),")",   &
                      "(",corn_lon(idx1,idx2+1),corn_lat(idx1,idx2+1),")"
            endif
          enddo
        enddo
        call flush(6)
      endif

      if ( verbose > 10 ) then
        !--- Print the grid mask (0 = cell not used, 1 = cell used)
        do idx1=ny,1,-1
          write(6,'(i4,2x,360i1)')idx1,local_imask(:,idx1)
        enddo
        call flush(6)
      endif

      !--- Initialize the grid_id cpl_grid element
      fidmax = 0
      if ( n_grid_list > 0 ) then
        !--- Determine the maximum file ID assigned to this grid description file
        do idx1=1,n_grid_list
          !--- Ignore grids that do not have file_name defined
          if ( len_trim(cpl_grid(idx1)%file_name) < 1 ) cycle
          !--- Extract the file name and file ID from the current cpl_grid file name
          strng = grid_desc_file_name(idx1, fid=fid)
          !--- If the cpl_grid file name matches the current file name then assign fidmax
          if ( trim(adjustl(strng)) .eq. trim(adjustl(file_name)) ) fidmax = max(fid,fidmax)
        enddo
      endif
      !--- The file ID for the current grid is one greater than the maximum file ID
      !--- that was found in cpl_grid for this file name or 1 if no grids matched
      fid = fidmax + 1
      sfx = " "
      write(sfx,'("(",i3.3,")")')fid

      cpl_grid(grid_id)%file_name = trim(file_name)//trim(sfx)
      cpl_grid(grid_id)%title     = trim(grid_title)
      cpl_grid(grid_id)%nx        = local_grid_dims(1)
      cpl_grid(grid_id)%ny        = local_grid_dims(2)
      cpl_grid(grid_id)%size      = local_grid_size
      cpl_grid(grid_id)%edges     = local_grid_edges

      write(6,'(2a)')"define_grid_from_file: file_name = ",trim(file_name)
      call flush(6)

      if ( cart ) then
        !--- Create a cartesian grid with no periodic dimensions
        ! 1st coord is 1D and depends on 1st Grid dim
        ! 2nd coord is 1D and depends on 2nd Grid dim
        grid_list(grid_id) = ESMF_GridCreateNoPeriDim(regDecomp=regDecomp,         &
                                                      maxIndex=(/nx,ny/),          &
                                                      coordSys=ESMF_COORDSYS_CART, &
                                                      coordDep1=(/one/),           &
                                                      coordDep2=(/two/),           &
                                                      indexflag=ESMF_INDEX_GLOBAL, &
                                                      name=trim(file_name),        &
                                                      rc=rc )
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Cannot create grid from "//trim(file_name), line=lnum, method=method)

      else
        !--- Create a grid containing 1 periodic dimension and a regular distribution
        !??? gridEdgeLWidth = (/ local_grid_edges(1,1), local_grid_edges(1,2) /)
        !??? gridEdgeUWidth = (/ local_grid_edges(2,1), local_grid_edges(2,2) /)
        gridEdgeLWidth = (/ 0, 0 /)
        gridEdgeUWidth = (/ 0, 1 /)
        grid_list(grid_id) = ESMF_GridCreate1PeriDim( regDecomp=regDecomp, &
                                                      maxIndex=(/nx,ny/), &
                                                      coordSys=ESMF_COORDSYS_SPH_DEG, &
                                                      gridEdgeLWidth=gridEdgeLWidth,  &
                                                      gridEdgeUWidth=gridEdgeUWidth,  &
                                                      indexflag=ESMF_INDEX_GLOBAL,    &
                                                      name=trim(file_name), &
                                                      rc=rc )
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Cannot create grid from "//trim(file_name), line=lnum, method=method)
      endif

      !--- Allocate storage for the grid center coordinates
      call ESMF_GridAddCoord(grid_list(grid_id), staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot allocate space for grid center coords", line=lnum, method=method)

      !--- Allocate storage for the grid corner coordinates
      call ESMF_GridAddCoord(grid_list(grid_id), staggerloc=ESMF_STAGGERLOC_CORNER, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot allocate space for grid corner coords", line=lnum, method=method)

      !--- Allocate storage for the grid mask
      call ESMF_GridAddItem(grid_list(grid_id), itemflag=ESMF_GRIDITEM_MASK, &
                                            staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot allocate space for grid mask", line=lnum, method=method)

      !--- Retrieve a pointer to the grid mask
      !--- add assign these values using data read from the file
      call ESMF_GridGetItem(grid_list(grid_id), &
                            staggerloc=ESMF_STAGGERLOC_CENTER, &
                            itemflag=ESMF_GRIDITEM_MASK, &
                            farrayPtr=grid_mask, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot retrieve grid mask", line=lnum, method=method)

      nxc = size(grid_mask,dim=1)
      nyc = size(grid_mask,dim=2)
      grid_mask = local_imask(1:nxc,1:nyc)

      !--- Allocate and assign space for this mask in cpl_grid
      if ( associated(cpl_grid(grid_id)%mask) ) deallocate(cpl_grid(grid_id)%mask)
      allocate( cpl_grid(grid_id)%mask(nxc,nyc) )
      cpl_grid(grid_id)%mask = local_imask(1:nxc,1:nyc)

      if ( cart ) then
        !--- This is a cartesian grid with no periodic dimensions
        !--- 1st coord (lon) is 1D and depends on 1st Grid dimension
        !--- 2nd coord (lat) is 1D and depends on 2nd Grid dimension
        !--- Retrieve a pointer to the cell center lon coordinate values
        !--- add assign these values using data read from the file
        call ESMF_GridGetCoord(grid_list(grid_id), coordDim=one, localDE=zero, &
                               staggerloc=ESMF_STAGGERLOC_CENTER, &
                               farrayPtr=grid_mid_lon_1d, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve grid_mid_lon coord values", line=lnum, method=method)

        grid_mid_lon_1d = mid_lon_1d(1:size(grid_mid_lon_1d))

        !--- Retrieve a pointer to the cell center lat coordinate values
        !--- add assign these values using data read from the file
        call ESMF_GridGetCoord(grid_list(grid_id), coordDim=two, localDE=zero, &
                               staggerloc=ESMF_STAGGERLOC_CENTER, &
                               farrayPtr=grid_mid_lat_1d, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve grid_mid_lat coord values", line=lnum, method=method)

        grid_mid_lat_1d = mid_lat_1d(1:size(grid_mid_lat_1d))

        !--- Retrieve a pointer to the cell corner lon coordinate values
        !--- add assign these values using data read from the file
        call ESMF_GridGetCoord(grid_list(grid_id), coordDim=one, localDE=zero, &
                               staggerloc=ESMF_STAGGERLOC_CORNER, &
                               farrayPtr=grid_corn_lon_1d, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve grid_corn_lon coord values", line=lnum, method=method)

        grid_corn_lon_1d = corn_lon_1d(1:size(grid_corn_lon_1d))

        !--- Retrieve a pointer to the cell corner lat coordinate values
        !--- add assign these values using data read from the file
        call ESMF_GridGetCoord(grid_list(grid_id), coordDim=two, localDE=zero, &
                               staggerloc=ESMF_STAGGERLOC_CORNER, &
                               farrayPtr=grid_corn_lat_1d, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve grid_corn_lat coord values", line=lnum, method=method)

        grid_corn_lat_1d = corn_lat_1d(1:size(grid_corn_lat_1d))

      else
        !--- This is a spherical grid containing 1 periodic dimension (lon)
        !--- 1st coord (lon) is 2D and depends on both Grid dimensions
        !--- 2nd coord (lat) is 2D and depends on both Grid dimensions

        !--- Retrieve a pointer to the cell center lon (dim 1) coordinate values
        !--- add assign these values using data read from the file
        call ESMF_GridGetCoord(grid_list(grid_id), coordDim=one, localDE=zero, &
                               staggerloc=ESMF_STAGGERLOC_CENTER, &
                               farrayPtr=grid_mid_lon, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve grid_mid_lon coord values", line=lnum, method=method)

        nxc = size(grid_mid_lon,dim=1)
        nyc = size(grid_mid_lon,dim=2)
        grid_mid_lon = mid_lon(1:nxc,1:nyc)

        !--- Retrieve a pointer to the cell center lat (dim 2) coordinate values
        !--- add assign these values using data read from the file
        call ESMF_GridGetCoord(grid_list(grid_id), coordDim=two, localDE=zero, &
                               staggerloc=ESMF_STAGGERLOC_CENTER, &
                               farrayPtr=grid_mid_lat, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve grid_mid_lat coord values", line=lnum, method=method)

        nxc = size(grid_mid_lat,dim=1)
        nyc = size(grid_mid_lat,dim=2)
        grid_mid_lat = mid_lat(1:nxc,1:nyc)

        !--- Retrieve a pointer to the cell corner lon (dim 1) coordinate values
        !--- add assign these values using data read from the file
        call ESMF_GridGetCoord(grid_list(grid_id), coordDim=one, localDE=zero, &
                               staggerloc=ESMF_STAGGERLOC_CORNER, &
                               farrayPtr=grid_corn_lon, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve grid_corn_lon coord values", line=lnum, method=method)

        nxc = size(grid_corn_lon,dim=1)
        nyc = size(grid_corn_lon,dim=2)
        grid_corn_lon = corn_lon(1:nxc,1:nyc)

        !--- Retrieve a pointer to the cell corner lat (dim 2) coordinate values
        !--- add assign these values using data read from the file
        call ESMF_GridGetCoord(grid_list(grid_id), coordDim=two, localDE=zero, &
                               staggerloc=ESMF_STAGGERLOC_CORNER, &
                               farrayPtr=grid_corn_lat, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve grid_corn_lat coord values", line=lnum, method=method)

        nxc = size(grid_corn_lat,dim=1)
        nyc = size(grid_corn_lat,dim=2)
        grid_corn_lat = corn_lat(1:nxc,1:nyc)

      endif

      if ( add_area_grid_item ) then
        !--- Define a real kind=8 field using the current grid
        area_field = ESMF_FieldCreate(grid_list(grid_id), ESMF_TYPEKIND_R8, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Cannot create area field", line=lnum, method=method)

        !--- Compute src grid cell areas and store them in the area_field object
        call ESMF_FieldRegridGetArea(area_field, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Cannot compute area field", line=lnum, method=method)

        !--- Get a pointer to the 2D area array now in area_field
        if ( associated(area_data1) ) nullify( area_data1 )
        call ESMF_FieldGet( area_field, farrayPtr=area_data1, rc=rc )
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Cannot extract area field", line=lnum, method=method)

        if ( .false. ) then
          !--- Read the grid area from a file
          !--- If a grid area file is not available then the area calculated above is used
          if ( associated(area_data3) ) deallocate( area_data3 )
          if ( grid_id == id_ocn_grid ) then
            if ( len_trim(nemo_mesh_mask_file) > 0 ) then
              write(6,*)"define_grid_from_file: Reading grid area from ", &
                        trim(nemo_mesh_mask_file)
              call flush(6)
              call read_NEMO_mesh_mask_file(trim(nemo_mesh_mask_file), area=area_data3, mask=mask_data)
              !--- The grid cell area read from this file is in meters squared
              !--- Normalize to the unit sphere assuming a radus of 6371220 meters
              !--- Note: NEMO uses earth radius = 6371229 m whereas the AGCM uses 6371220 m
              !---       We want to use a common radius for both so we choose 6371220 m
              rad2 = 6371220.0_8 * 6371220.0_8
              area_data3 = area_data3/rad2
              garea = sum(area_data3)
              write(6,*)"define_grid_from_file: Normalized NEMO total grid area ",garea
              write(6,*)"define_grid_from_file: NEMO mask shape ",shape(mask_data), &
                        "  minval=",minval(mask_data),"  maxval=",maxval(mask_data), &
                        "  num==0 ",count(mask_data == 0.0_8),"  num==1 ",count(mask_data == 1.0_8)
            else
              write(6,*)"define_grid_from_file: NEMO grid area file is not available."
              call flush(6)
            endif
          else if ( grid_id == id_atm_grid ) then
            if ( len_trim(agcm_grid_area_file) > 0 ) then
              write(6,*)"define_grid_from_file: Reading grid area from ", &
                        trim(agcm_grid_area_file)
              call flush(6)
              call read_AGCM_grid_area_file(trim(agcm_grid_area_file), area=area_data3)
              !--- The grid cell area read from this file is in meters squared
              !--- Normalize to the unit sphere assuming a radus of 6371220 meters
              !--- Note: NEMO uses earth radius = 6371229 m whereas the AGCM uses 6371220 m
              !---       We want to use a common radius for both so we choose 6371220 m
              rad2 = 6371220.0_8 * 6371220.0_8
              area_data3 = area_data3/rad2
              garea = sum(area_data3)
              write(6,*)"define_grid_from_file: Normalized AGCM total grid area ",garea
            else
              write(6,*)"define_grid_from_file: AGCM grid area file is not available."
              call flush(6)
            endif
          endif
          if ( associated(area_data3) ) then
            if ( size(area_data3,dim=1) /= size(area_data1,dim=1) .or. &
                 size(area_data3,dim=2) /= size(area_data1,dim=2) ) then
              write(6,*)"define_grid_from_file: Shape mismatch for grid area read from file"
              write(6,*)"   found nx,ny = ",size(area_data3,dim=1),size(area_data3,dim=2)
              write(6,*)"expected nx,ny = ",size(area_data1,dim=1),size(area_data1,dim=2)
              call xit("define_grid_from_file",-1)
            endif
            write(6,'(a,g24.16)')"define_grid_from_file: ESMF grid total area is ",sum(area_data1)
            write(6,'(a,g24.16)')"define_grid_from_file: File grid total area is ",sum(area_data3)
            call flush(6)
            !--- Overwrite area calculated via ESMF_FieldRegridGetArea with area read from a file
            area_data1 = area_data3
            deallocate(area_data3)
          endif
        endif

        if ( grid_id == id_atm_grid ) then
          !--- Assign this cell area data to the global array atm_grid_cell_area
          if ( associated(atm_grid_cell_area) ) deallocate(atm_grid_cell_area)
          allocate( atm_grid_cell_area(size(area_data1,dim=1), size(area_data1,dim=2)) )
          atm_grid_cell_area = area_data1
        else if ( grid_id == id_ocn_grid ) then
          !--- Assign this cell area data to the global array ocn_grid_cell_area
          if ( associated(ocn_grid_cell_area) ) deallocate(ocn_grid_cell_area)
          allocate( ocn_grid_cell_area(size(area_data1,dim=1), size(area_data1,dim=2)) )
          ocn_grid_cell_area = area_data1
        endif

        !--- Allocate space for the area grid item in this grid
        call ESMF_GridAddItem(grid_list(grid_id), itemflag=ESMF_GRIDITEM_AREA, &
             staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Cannot allocate space for grid area", line=lnum, method=method)

        !--- Retrieve a pointer to the grid area in the current grid
        call ESMF_GridGetItem(grid_list(grid_id), &
                              staggerloc=ESMF_STAGGERLOC_CENTER, &
                              itemflag=ESMF_GRIDITEM_AREA, &
                              farrayPtr=area_data2, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Cannot retrieve grid area", line=lnum, method=method)

        !--- Verify that the referenced array sizes are equal
        if ( size(area_data1,dim=1) /= size(area_data2,dim=1) .or. &
             size(area_data1,dim=2) /= size(area_data2,dim=2) ) then
          write(6,*)"define_grid_from_file: Shape mismatch for grid areas"
          write(6,*)"   found nx,ny = ",size(area_data1,dim=1),size(area_data1,dim=2)
          write(6,*)"expected nx,ny = ",size(area_data2,dim=1),size(area_data2,dim=2)
          call xit("define_grid_from_file",-1)
        endif

        !--- Assign the grid area in the current grid with the data from area_field
        area_data2 = area_data1

        if (verbose > -10) then
          write(6,'(a,g24.16)')"define_grid_from_file: grid total area is ",sum(area_data1)
          call flush(6)
        endif

        !--- Clean up area arrays
        nullify(area_data1)
        nullify(area_data2)
        call ESMF_FieldDestroy(area_field, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Problem destroying area_field", line=lnum, method=method)
      endif

!KEEP-FOR-NOW ESMF_GridCreate does almost exactly what this subroutine does with a few exceptions
!KEEP-FOR-NOW such as handling overlap longitudes in the SCRIP grid description file,
!KEEP-FOR-NOW specifying the grid name and optionally adding extra ghost cells
!KEEP-FOR-NOW
!KEEP-FOR-NOW      !--- Create an ESMF grid object for the default atm grid using data read
!KEEP-FOR-NOW      !--- from a SCRIP format netcdf file
!KEEP-FOR-NOW      grid_list(grid_id) = ESMF_GridCreate(filename=trim(file_name), &
!KEEP-FOR-NOW                                           fileFormat=ESMF_FILEFORMAT_SCRIP, &
!KEEP-FOR-NOW                                           regDecomp=regDecomp, &
!KEEP-FOR-NOW                                           isSphere=T_esmf, &
!KEEP-FOR-NOW                                           addCornerStagger=T_esmf, &
!KEEP-FOR-NOW                                           rc=rc)
!KEEP-FOR-NOW      call handle_ESMF_err(rc, "Cannot create grid from "//trim(file_name))

      !--- Associate this grid with the cpl_grid array
      cpl_grid(grid_id)%grid => grid_list(grid_id)

! Do not release this file because other routines will need to read it
! ??? Should it be an option to remove or not remove ???
!      !--- Clean up
!      if (file_was_accessed) then
!        !--- Release this file, it was aquired by access
!        write(6,'(2a)')"define_grid_from_file: Release remote file ",trim(file_name)
!        strng=' '
!        strng="release "//trim(file_name)
!        call system(trim(strng))
!      endif

      if ( associated(local_imask) )    deallocate(local_imask)
      if ( associated(mid_lon) )  deallocate(mid_lon)
      if ( associated(mid_lat) )  deallocate(mid_lat)
      if ( associated(corn_lon) ) deallocate(corn_lon)
      if ( associated(corn_lat) ) deallocate(corn_lat)
      if ( associated(mid_lon_1d) )  deallocate(mid_lon_1d)
      if ( associated(mid_lat_1d) )  deallocate(mid_lat_1d)
      if ( associated(corn_lon_1d) ) deallocate(corn_lon_1d)
      if ( associated(corn_lat_1d) ) deallocate(corn_lat_1d)

    end subroutine define_grid_from_file

    !***************************************************************************
    !--- Compute grid cell area using ESMF_FieldRegridGetArea
    !--- Given a grid ID (ie index into the global grid_list array) compute
    !--- grid cell area for the grid defined at that index
    !***************************************************************************
    subroutine compute_grid_area(gridID, area)
      !--- The index in grid_list of the grid for which the area is requested
      integer, intent(in) :: gridID

      real(kind=8), pointer :: area(:,:)

      !--- Local
      character(64) :: method = "compute_grid_area"
      integer(ESMF_KIND_I) :: lnum
      integer(ESMF_KIND_I) :: rc
      type(ESMF_Field) :: area_field
      real(ESMF_KIND_R8), pointer :: area_data(:,:) => null()

      !--- Calculate the grid cell area internally
      !--- Define a real kind=8 field using the current grid
      area_field = ESMF_FieldCreate(grid_list(gridID), ESMF_TYPEKIND_R8, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot create area field", line=lnum, method=method)

      !--- Compute src grid cell areas and store them in the area_field object
      call ESMF_FieldRegridGetArea(area_field, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot compute area field", line=lnum, method=method)

      !--- Get a pointer to the 2D area array now in area_field
      if ( associated(area_data) ) nullify( area_data )
      call ESMF_FieldGet( area_field, farrayPtr=area_data, rc=rc )
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot extract area field", line=lnum, method=method)

      if ( associated(area) ) deallocate(area)
      allocate( area(size(area_data,dim=1),size(area_data,dim=2)) )
      area = area_data

      call ESMF_FieldDestroy(area_field, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Problem destroying area_field", line=lnum, method=method)

    end subroutine compute_grid_area

    !***************************************************************************
    !--- Determine if all relevant AGCM grid information is available in memory.
    !--- If the AGCM is coupled, this infomation is passed from the AGCM to
    !--- the coupler at runtime. If it is not then this information is read from
    !--- a file named cplrs_atm_grid.nc which is saved in the coupler restart.
    !--- The following required variables are found in the module cpl_com
    !---    atm_lon, atm_lat, atm_gwt, atm_grid_cell_area, atm_nlon, atm_nlat
    !---    atm_earth_radius
    !***************************************************************************
    function agcm_coords_available() result(valid)
      logical :: valid

      valid = .true.

      !--- If any variables are undefined then reset valid as false
      if ( .not. associated(atm_lon) .or. &
           .not. associated(atm_lat) .or. &
           .not. associated(atm_gwt) .or. &
           .not. associated(atm_grid_cell_area)  ) then
        write(6,*)"agcm_coords_available: Incomplete AGCM grid info:"
        if ( .not. associated(atm_lon) ) write(6,*)" atm_lon is missing"
        if ( .not. associated(atm_lat) ) write(6,*)" atm_lat is missing"
        if ( .not. associated(atm_gwt) ) write(6,*)" atm_gwt is missing"
        if ( .not. associated(atm_grid_cell_area) ) write(6,*)" atm_grid_cell_area is missing"
        valid = .false.
        return
      endif

      if ( atm_nlon < 0 ) then
        write(6,*)"agcm_coords_available: atm_nlon is not defined"
        valid = .false.
      endif

      if ( atm_nlat < 0 ) then
        write(6,*)"agcm_coords_available: atm_nlat is not defined"
        valid = .false.
      endif

      if ( atm_earth_radius < 0.0 ) then
        write(6,*)"agcm_coords_available: atm_earth_radius is not defined"
        valid = .false.
      endif

      if ( size(atm_lon,dim=1) /= atm_nlon ) then
        write(6,*)"agcm_coords_available: Size mismatch for AGCM lon coord"
        write(6,*)"   expected: ",atm_nlon
        write(6,*)"      found: ",shape(atm_lon)
        valid = .false.
      endif

      if ( size(atm_lat,dim=1) /= atm_nlat ) then
        write(6,*)"agcm_coords_available: Size mismatch for AGCM lat coord"
        write(6,*)"   expected: ",atm_nlat
        write(6,*)"      found: ",shape(atm_lat)
        valid = .false.
      endif

      if ( size(atm_gwt,dim=1) /= atm_nlat ) then
        write(6,*)"agcm_coords_available: Size mismatch for AGCM Gaussian weight"
        write(6,*)"   expected: ",atm_nlat
        write(6,*)"      found: ",shape(atm_gwt)
        valid = .false.
      endif

      if ( size(atm_grid_cell_area,dim=1) /= atm_nlon .or. &
           size(atm_grid_cell_area,dim=2) /= atm_nlat ) then
        write(6,*)"agcm_coords_available: Shape mismatch for AGCM grid cell area"
        write(6,*)"   expected: ",atm_nlon,atm_nlat
        write(6,*)"      found: ",shape(atm_grid_cell_area)
        valid = .false.
      endif

    end function agcm_coords_available

    !***************************************************************************
    !--- Define a grid object using internal AGCM coordinate and mask data
    !--- This information is passed to the coupler from the AGCM at run time
    !***************************************************************************
    subroutine define_grid_from_agcm_coords(grid_id)
      !--- The index in grid_list and cpl_grid of this grid
      integer, intent(in) :: grid_id

      !--- Local
      character(64) :: method = "define_grid_from_agcm_coords"
      integer(ESMF_KIND_I) :: lnum
      integer(ESMF_KIND_I) :: rc
      integer(ESMF_KIND_I) :: regDecomp(2)
      integer(ESMF_KIND_I) :: nx_esmf, ny_esmf
      integer(ESMF_KIND_I) :: gridEdgeLWidth(2), gridEdgeUWidth(2)
      type(ESMF_DistGrid) :: distGrid
      logical :: exists, use_atm_grid_cell_area
      integer :: idx1, idx2, fid, fidmax, xdim(6), ydim(6)
      integer :: nx, ny, ix, iy
      integer :: nxc, nyc
      integer :: offset(2)
      character(128) :: sfx, strng
      real(kind=8)   :: gsize
      integer(kind=4) :: ncfid

      character(128)         :: grid_title
      real(ESMF_KIND_R8), pointer, dimension(:,:) :: mid_lon=>null(), mid_lat=>null()
      real(ESMF_KIND_R8), pointer, dimension(:,:) :: corn_lon=>null(), corn_lat=>null()
      real(ESMF_KIND_R8), pointer, dimension(:,:) :: cell_area=>null()
      real(ESMF_KIND_R8), pointer :: lonm(:) => null()
      real(ESMF_KIND_R8), pointer :: lonb(:) => null()
      real(ESMF_KIND_R8), pointer :: latm(:) => null()
      real(ESMF_KIND_R8), pointer :: latb(:) => null()
      real(ESMF_KIND_R8) :: delta

      !--- Internal info associated with this grid
      real(ESMF_KIND_R8), pointer    :: grid_mid_lon(:,:) => null()
      real(ESMF_KIND_R8), pointer    :: grid_mid_lat(:,:) => null()
      real(ESMF_KIND_R8), pointer    :: grid_corn_lon(:,:) => null()
      real(ESMF_KIND_R8), pointer    :: grid_corn_lat(:,:) => null()
      integer(ESMF_KIND_I4), pointer :: grid_mask(:,:) => null()

      type(ESMF_Field) :: area_field
      real(ESMF_KIND_R8), pointer :: area_data(:,:) => null()
      real(ESMF_KIND_R8), pointer :: mask_data(:,:) => null()
      real(ESMF_KIND_R8) :: garea, rad2, sum_gwt

      if ( .not. agcm_coords_available() ) then
        !--- AGCM coords are not defined, likely because the AGCM is not coupled
        !--- In this case populate AGCM grid infomation arrays with data read from a file
        !--- This file is part of the coupler restart and so should be available at run time
        call read_atm_grid_file("cplrs_atm_grid.nc")
        if ( .not. agcm_coords_available() ) then
          write(6,*)"define_grid_from_agcm_coords: AGCM coordinate information is not available"
          call xit("define_grid_from_agcm_coords",-1)
        endif
      endif

      !--- Set a grid title for use with cpl_grid(grid_id)
      grid_title = "AGCM Grid"

      !--- Define 1D versions of lon/lat cell centers and boundaries
      allocate( lonm(atm_nlon), lonb(atm_nlon+1) )
      allocate( latm(atm_nlat), latb(atm_nlat+1) )

      !--- Values in atm_lon represent grid cell centers
      lonm = atm_lon

      !--- Longitude is equally spaced so boundaries are center-delta, center+delta
      !--- where delta is one half the constant grid cell width
      delta = 180.0_8/real(atm_nlon,kind=8)
      do idx1=1,atm_nlon
        lonb(idx1) = atm_lon(idx1) - delta
      enddo
      lonb(atm_nlon+1) = atm_lon(atm_nlon) + delta

      !--- Values in atm_lat represent grid cell centers
      latm = atm_lat

      !--- Calculate boundaries of gaussian latitude cell centers
      !--- End points are the North and South poles
      !--- The equator is always a boundary
      latb(1) = -90.0_8
      latb(atm_nlat/2+1) = 0.0_8
      sum_gwt = -1.0_8
      do idx1=2,atm_nlat/2
        sum_gwt = sum_gwt + atm_gwt(idx1-1)
        latb(idx1) = asin(sum_gwt) * rad2deg
        latb(atm_nlat-idx1+2) = -latb(idx1)
      enddo
      !--- Assign the North pole value to be exactly 90
      latb(atm_nlat+1) = 90.0_8

      !--- Cell center longitude
      if ( associated(mid_lon) ) deallocate(mid_lon)
      allocate( mid_lon(atm_nlon,atm_nlat) )
      !--- Broadcast the 1D array lonm accross the 2D mid_lon array
      do idx1=1,atm_nlat
        mid_lon(:,idx1) = lonm(1:atm_nlon)
      enddo

      !--- Cell center latitude
      if ( associated(mid_lat) ) deallocate(mid_lat)
      allocate( mid_lat(atm_nlon,atm_nlat) )
      !--- Broadcast the 1D array latm accross the 2D mid_lat array
      do idx1=1,atm_nlon
        mid_lat(idx1,:) = latm(1:atm_nlat)
      enddo

      !--- Cell corner longitude
      if ( associated(corn_lon) ) deallocate(corn_lon)
      allocate( corn_lon(atm_nlon+1,atm_nlat+1) )
      !--- Broadcast the 1D array lonb accross the 2D corn_lon array
      do idx1=1,atm_nlat+1
        corn_lon(:,idx1) = lonb(1:atm_nlon+1)
      enddo

      !--- Cell corner latitude
      if ( associated(corn_lat) ) deallocate(corn_lat)
      allocate( corn_lat(atm_nlon+1,atm_nlat+1) )
      !--- Broadcast the 1D array latb accross the 2D corn_lat array
      do idx1=1,atm_nlon+1
        corn_lat(idx1,:) = latb(1:atm_nlat+1)
      enddo

      !--- Define regular decomposition counts for each dimension
      regDecomp(1) = cpl_PetCount
      regDecomp(2) = 1

      !--- Initialize the grid_id cpl_grid element
      fidmax = 0
      if ( n_grid_list > 0 ) then
        !--- Determine the maximum file ID assigned to this grid description file
        do idx1=1,n_grid_list
          !--- Ignore grids that do not have file_name defined
          if ( len_trim(cpl_grid(idx1)%file_name) < 1 ) cycle
          !--- Extract the file name and file ID from the current cpl_grid file name
          strng = grid_desc_file_name(idx1, fid=fid)
          !--- If the cpl_grid file name matches the current file name then assign fidmax
          if ( trim(adjustl(strng)) .eq. "internal" ) fidmax = max(fid,fidmax)
        enddo
      endif
      !--- The file ID for the current grid is one greater than the maximum file ID
      !--- that was found in cpl_grid for this file name or 1 if no grids matched
      fid = fidmax + 1
      sfx = " "
      write(sfx,'("(",i3.3,")")')fid

      cpl_grid(grid_id)%file_name = "internal"//trim(sfx)
      cpl_grid(grid_id)%title     = trim(grid_title)
      cpl_grid(grid_id)%nx        = atm_nlon
      cpl_grid(grid_id)%ny        = atm_nlat
      cpl_grid(grid_id)%size      = atm_nlon*atm_nlat
      cpl_grid(grid_id)%edges     = 0

      !--- Create a grid containing 1 periodic dimension and a regular distribution
      gridEdgeLWidth = (/ 0, 0 /)
      gridEdgeUWidth = (/ 0, 1 /)
      nx_esmf = atm_nlon
      ny_esmf = atm_nlat
      grid_list(grid_id) = ESMF_GridCreate1PeriDim( regDecomp=regDecomp, &
                                                    maxIndex=(/nx_esmf,ny_esmf/), &
                                                    coordSys=ESMF_COORDSYS_SPH_DEG, &
                                                    gridEdgeLWidth=gridEdgeLWidth,  &
                                                    gridEdgeUWidth=gridEdgeUWidth,  &
                                                    indexflag=ESMF_INDEX_GLOBAL,    &
                                                    name=trim(grid_title), &
                                                    rc=rc )
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot create grid from NEMO coordinates", line=lnum, method=method)

      !--- Allocate storage for the grid center coordinates
      call ESMF_GridAddCoord(grid_list(grid_id), staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot allocate space for grid center coords", line=lnum, method=method)

      !--- Allocate storage for the grid corner coordinates
      call ESMF_GridAddCoord(grid_list(grid_id), staggerloc=ESMF_STAGGERLOC_CORNER, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot allocate space for grid corner coords", line=lnum, method=method)

      !--- Allocate storage for the grid mask
      call ESMF_GridAddItem(grid_list(grid_id), itemflag=ESMF_GRIDITEM_MASK, &
                                            staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot allocate space for grid mask", line=lnum, method=method)

      !--- Retrieve a pointer to the grid mask
      !--- add assign these values using data read from the file
      call ESMF_GridGetItem(grid_list(grid_id), &
                            staggerloc=ESMF_STAGGERLOC_CENTER, &
                            itemflag=ESMF_GRIDITEM_MASK, &
                            farrayPtr=grid_mask, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot retrieve grid mask", line=lnum, method=method)

      nxc = size(grid_mask,dim=1)
      nyc = size(grid_mask,dim=2)
      !--- The AGCM integer mask is always 1 everywhere
      grid_mask = 1

      !--- Allocate and assign space for this mask in cpl_grid
      if ( associated(cpl_grid(grid_id)%mask) ) deallocate(cpl_grid(grid_id)%mask)
      allocate( cpl_grid(grid_id)%mask(nxc,nyc) )
      !--- The AGCM integer mask is always 1 everywhere
      cpl_grid(grid_id)%mask = 1

      !--- This is a spherical grid containing 1 periodic dimension (lon)
      !--- 1st coord (lon) is 2D and depends on both Grid dimensions
      !--- 2nd coord (lat) is 2D and depends on both Grid dimensions

      !--- Retrieve a pointer to the cell center lon (dim 1) coordinate values
      !--- add assign these values using data read from the file
      call ESMF_GridGetCoord(grid_list(grid_id), coordDim=one, localDE=zero, &
                             staggerloc=ESMF_STAGGERLOC_CENTER, &
                             farrayPtr=grid_mid_lon, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to retrieve grid_mid_lon coord values", line=lnum, method=method)

      nxc = size(grid_mid_lon,dim=1)
      nyc = size(grid_mid_lon,dim=2)
      grid_mid_lon = mid_lon(1:nxc,1:nyc)

      !--- Retrieve a pointer to the cell center lat (dim 2) coordinate values
      !--- add assign these values using data read from the file
      call ESMF_GridGetCoord(grid_list(grid_id), coordDim=two, localDE=zero, &
                             staggerloc=ESMF_STAGGERLOC_CENTER, &
                             farrayPtr=grid_mid_lat, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to retrieve grid_mid_lat coord values", line=lnum, method=method)

      nxc = size(grid_mid_lat,dim=1)
      nyc = size(grid_mid_lat,dim=2)
      grid_mid_lat = mid_lat(1:nxc,1:nyc)

      !--- Retrieve a pointer to the cell corner lon (dim 1) coordinate values
      !--- add assign these values using data read from the file
      call ESMF_GridGetCoord(grid_list(grid_id), coordDim=one, localDE=zero, &
                             staggerloc=ESMF_STAGGERLOC_CORNER, &
                             farrayPtr=grid_corn_lon, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to retrieve grid_corn_lon coord values", line=lnum, method=method)

      nxc = size(grid_corn_lon,dim=1)
      nyc = size(grid_corn_lon,dim=2)
      grid_corn_lon = corn_lon(1:nxc,1:nyc)

      !--- Retrieve a pointer to the cell corner lat (dim 2) coordinate values
      !--- add assign these values using data read from the file
      call ESMF_GridGetCoord(grid_list(grid_id), coordDim=two, localDE=zero, &
                             staggerloc=ESMF_STAGGERLOC_CORNER, &
                             farrayPtr=grid_corn_lat, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to retrieve grid_corn_lat coord values", line=lnum, method=method)

      nxc = size(grid_corn_lat,dim=1)
      nyc = size(grid_corn_lat,dim=2)
      grid_corn_lat = corn_lat(1:nxc,1:nyc)

      !--- Grid cell area
      write(6,*)"define_grid_from_agcm_coords: Dimensional AGCM total grid area ",sum(atm_grid_cell_area)
      use_atm_grid_cell_area = .true.
      if ( use_atm_grid_cell_area ) then
        !--- Use the grid cell area defined in the AGCM and passed to the coupler
        if ( associated(cell_area) ) deallocate(cell_area)
        allocate( cell_area(atm_nlon,atm_nlat) )
        cell_area = atm_grid_cell_area

        !--- Renormalize the area to the unit sphere for use in regridding
        rad2 = atm_earth_radius * atm_earth_radius
        cell_area = cell_area/rad2
        write(6,*)"define_grid_from_agcm_coords: Normalized total ocean area on AGCM grid ",sum(cell_area)

      else
        !--- Calculate grid cell area internally
        !--- Define a real kind=8 field using the current grid
        call compute_grid_area(grid_id, area_data)

        if ( atm_nlon /= size(area_data,dim=1) .or. &
             atm_nlat /= size(area_data,dim=2) ) then
          write(6,*)"define_grid_from_agcm_coords: Shape mismatch for grid areas"
          write(6,*)"expected nx,ny = ",atm_nlon,atm_nlat
          write(6,*)"   found nx,ny = ",size(area_data,dim=1),size(area_data,dim=2)
          call xit("define_grid_from_agcm_coords",-1)
        endif

        !--- Assign this calculated area to the cell_area array
        if ( associated(cell_area) ) deallocate(cell_area)
        allocate( cell_area(atm_nlon,atm_nlat) )
        cell_area = area_data
      endif

      !--- Allocate space for the area grid item in this grid
      call ESMF_GridAddItem(grid_list(grid_id), itemflag=ESMF_GRIDITEM_AREA, &
           staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot allocate space for grid area", line=lnum, method=method)

      !--- Retrieve a pointer to the grid area in the current grid
      nullify(area_data)
      call ESMF_GridGetItem(grid_list(grid_id), &
                            staggerloc=ESMF_STAGGERLOC_CENTER, &
                            itemflag=ESMF_GRIDITEM_AREA, &
                            farrayPtr=area_data, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot retrieve grid area", line=lnum, method=method)

      !--- Verify that the referenced array sizes are equal
      if ( size(cell_area,dim=1) /= size(area_data,dim=1) .or. &
           size(cell_area,dim=2) /= size(area_data,dim=2) ) then
        write(6,*)"define_grid_from_agcm_coords: Shape mismatch for grid areas"
        write(6,*)"expected nx,ny = ",size(cell_area,dim=1),size(cell_area,dim=2)
        write(6,*)"   found nx,ny = ",size(area_data,dim=1),size(area_data,dim=2)
        call xit("define_grid_from_agcm_coords",-1)
      endif

      !--- Assign the grid area in the current grid with the data from cell_area
      area_data = cell_area

      if (verbose > -10) then
        write(6,*)"define_grid_from_agcm_coords:    Internal AGCM total grid area ",sum(area_data)
        call flush(6)
      endif

      !--- Associate this grid with the cpl_grid array
      cpl_grid(grid_id)%grid => grid_list(grid_id)

      !--- Clean up area arrays
      nullify(area_data)
      if ( associated(mid_lon) )  deallocate(mid_lon)
      if ( associated(mid_lat) )  deallocate(mid_lat)
      if ( associated(corn_lon) ) deallocate(corn_lon)
      if ( associated(corn_lat) ) deallocate(corn_lat)
      if ( associated(cell_area)) deallocate(cell_area)
      if ( associated(lonm)) deallocate(lonm)
      if ( associated(lonb)) deallocate(lonb)
      if ( associated(latm)) deallocate(latm)
      if ( associated(latb)) deallocate(latb)

    end subroutine define_grid_from_agcm_coords

    !***************************************************************************
    !--- Determine if all relevant NEMO grid information is available in memory.
    !--- If NEMO is coupled, this infomation is passed from NEMO to the coupler
    !--- at runtime. If it is not then this information is read from a file
    !--- named cplrs_ocn_grid.nc which is saved in the coupler restart.
    !--- The following required variables are found in the module cpl_com
    !---    nemo_glamt, nemo_gphit, nemo_glamf, nemo_gphif, nemo_tmask, 
    !---    nemo_e1t, nemo_e2t
    !***************************************************************************
    function nemo_coords_available() result(valid)
      logical :: valid

      !--- Local
      integer :: idx1, nx, ny, xdim(6), ydim(6)

      valid = .true.

      !--- If any required grid information is missing then set valid false
      if ( .not. associated(nemo_glamt) .or. &
           .not. associated(nemo_gphit) .or. &
           .not. associated(nemo_glamf) .or. &
           .not. associated(nemo_gphif) .or. &
           .not. associated(nemo_tmask) .or. &
           .not. associated(nemo_e1t)   .or. &
           .not. associated(nemo_e2t)  ) then
        write(6,*)"nemo_coords_available: NEMO coordinate info is not avaialble"
        if ( .not. associated(nemo_glamt) ) write(6,*)"    nemo_glamt is missing"
        if ( .not. associated(nemo_gphit) ) write(6,*)"    nemo_gphit is missing"
        if ( .not. associated(nemo_glamf) ) write(6,*)"    nemo_glamf is missing"
        if ( .not. associated(nemo_gphif) ) write(6,*)"    nemo_gphif is missing"
        if ( .not. associated(nemo_tmask) ) write(6,*)"    nemo_tmask is missing"
        if ( .not. associated(nemo_e1t)   ) write(6,*)"    nemo_e1t   is missing"
        if ( .not. associated(nemo_e2t)   ) write(6,*)"    nemo_e2t   is missing"
        valid = .false.
        return
      endif

      !--- Verify all variables have identical shape
      nx = size(nemo_glamt,dim=1)
      ny = size(nemo_glamt,dim=2)
      xdim(1) = size(nemo_gphit,dim=1)
      ydim(1) = size(nemo_gphit,dim=2)
      xdim(2) = size(nemo_glamf,dim=1)
      ydim(2) = size(nemo_glamf,dim=2)
      xdim(3) = size(nemo_gphif,dim=1)
      ydim(3) = size(nemo_gphif,dim=2)
      xdim(4) = size(nemo_tmask,dim=1)
      ydim(4) = size(nemo_tmask,dim=2)
      xdim(5) = size(nemo_e1t,dim=1)
      ydim(5) = size(nemo_e1t,dim=2)
      xdim(6) = size(nemo_e2t,dim=1)
      ydim(6) = size(nemo_e2t,dim=2)
      if ( xdim(1) /= nx .or. ydim(1) /= ny .or. &
           xdim(2) /= nx .or. ydim(2) /= ny .or. &
           xdim(3) /= nx .or. ydim(3) /= ny .or. &
           xdim(4) /= nx .or. ydim(4) /= ny .or. &
           xdim(5) /= nx .or. ydim(5) /= ny .or. &
           xdim(6) /= nx .or. ydim(6) /= ny) then
        write(6,*)"nemo_coords_available: Size mismatch for NEMO coordinates"
        write(6,*)"   expected shape: ",nx,ny
        do idx1=1,6
          write(6,*)idx1," found shape: ",xdim(idx1),ydim(idx1)
        enddo
        valid = .false.
      endif

    end function nemo_coords_available

    !***************************************************************************
    !--- Define a grid object using internal NEMO coordinate and mask data
    !--- This information is passed to the coupler from NEMO at run time
    !***************************************************************************
    subroutine define_grid_from_nemo_coords(grid_id)
      !--- The index in grid_list and cpl_grid of this grid
      integer, intent(in) :: grid_id

      !--- Local
      character(64) :: method = "define_grid_from_nemo_coords"
      integer(ESMF_KIND_I) :: lnum
      integer(ESMF_KIND_I) :: rc
      integer(ESMF_KIND_I) :: regDecomp(2)
      integer(ESMF_KIND_I) :: nx_esmf, ny_esmf
      integer(ESMF_KIND_I) :: gridEdgeLWidth(2), gridEdgeUWidth(2)
      type(ESMF_DistGrid) :: distGrid
      logical :: exists, use_ocn_grid_cell_area
      integer :: idx1, idx2, fid, fidmax
      integer :: nx, ny, ic, jc
      integer :: nxc, nyc
      integer :: offset(2)
      character(128) :: sfx, strng
      real(kind=8)   :: gsize
      integer(kind=4) :: ncfid

      !--- Info to be read from the grid description file
      character(1)           :: gtype
      character(128)         :: grid_title
      integer,      pointer, dimension(:,:) :: imask=>null()
      real(kind=8), pointer, dimension(:,:) :: mid_lon=>null(), mid_lat=>null()
      real(kind=8), pointer, dimension(:,:) :: corn_lon=>null(), corn_lat=>null()
      real(kind=8), pointer, dimension(:,:) :: clon=>null(), clat=>null()
      real(kind=8), pointer, dimension(:,:) :: e1=>null(), e2=>null(), area=>null()
      real(kind=8), pointer :: awrk2d(:,:)=>null(), bwrk2d(:,:)=>null()
      integer,      pointer :: iwrk2d(:,:)=>null()

      !--- Internal info associated with this grid
      real(ESMF_KIND_R8), pointer    :: grid_mid_lon(:,:) => null()
      real(ESMF_KIND_R8), pointer    :: grid_mid_lat(:,:) => null()
      real(ESMF_KIND_R8), pointer    :: grid_corn_lon(:,:) => null()
      real(ESMF_KIND_R8), pointer    :: grid_corn_lat(:,:) => null()
      integer(ESMF_KIND_I4), pointer :: grid_mask(:,:) => null()

      !--- Grid area related
      type(ESMF_Field) :: area_field
      real(ESMF_KIND_R8), pointer :: area_data1(:,:) => null()
      real(ESMF_KIND_R8), pointer :: area_data2(:,:) => null()
      real(ESMF_KIND_R8), pointer :: area_data3(:,:) => null()
      real(ESMF_KIND_R8), pointer :: mask_data(:,:) => null()
      real(ESMF_KIND_R8) :: garea, rad2

      if ( .not. nemo_coords_available() ) then
        !--- NEMO coords are not defined, likely because NEMO is not coupled for this run
        !--- In this case populate NEMO grid infomation arrays with data read from a file
        !--- This file is part of the coupler restart and so should be available at run time
        call read_ocn_grid_file("cplrs_ocn_grid.nc")
        if ( .not. nemo_coords_available() ) then
          write(6,*)"define_grid_from_nemo_coords: NEMO coordinate information is not available"
          call xit("define_grid_from_nemo_coords",-1)
        endif
      endif

      gtype = "T"

      !--- Define variables that correspond to grid cell centers, corners and scale factors
      !--- offset is what needs to be added to (i,j) to get to the correct elements
      !--- in the corner arrays that correspond to the northeast corner of each cell
      select case(gtype)
        case ("t", "T")
          !--- The T grid requires
          !--- nemo_glamt
          !--- nemo_gphit
          !--- nemo_glamf
          !--- nemo_gphif
          !--- nemo_tmask
          !--- nemo_e1t
          !--- nemo_e2t
          !--- offset = (/0,0/)
          gtype = 'T'

          !--- Set a grid title for use with cpl_grid(grid_id)
          grid_title = "ORCA1 T Grid"

          nx = size(nemo_glamt,dim=1)
          ny = size(nemo_glamt,dim=2)

          if ( associated(mid_lon) ) deallocate(mid_lon)
          allocate( mid_lon(nx,ny) )
          !--- Grid cell center longitude
          mid_lon = nemo_glamt
          !--- Remove overlap longitudes
          call resize_nemo_field(mid_lon)

          !--- Grid cell center latitude
          if ( associated(mid_lat) ) deallocate(mid_lat)
          allocate( mid_lat(nx,ny) )
          mid_lat = nemo_gphit
          !--- Remove overlap longitudes
          call resize_nemo_field(mid_lat)

          !--- Longitude of north west corner of grid cell centered at (mid_lon,mid_lat)
          if ( associated(clon) ) deallocate(clon)
          allocate( clon(nx,ny) )
          clon = nemo_glamf

          !--- Latitude of north west corner of grid cell centered at (mid_lon,mid_lat)
          if ( associated(clat) ) deallocate(clat)
          allocate( clat(nx,ny) )
          clat = nemo_gphif

          !--- Reformat these corner lon/lat values to add missing boundary points
          !--- corn_lon/corn_lat are returned with 1 more row and column than clon/clat
          !--- Upon return (corn_lon(i,j),corn_lat(i,j)) corresponds to the south west
          !--- corner of a box centered at (mid_lon(i,j),mid_lat(i,j))
          !--- Put another way, the 4 corners have indicies
          !--- SW corner at (i  ,j  )
          !--- SE corner at (i+1,j  )
          !--- NW corner at (i  ,j+1)
          !--- NE corner at (i+1,j+1)
          offset = (/0,0/)
          call mm_corner_lon(clon, offset, corn_lon)
          call mm_corner_lat(clat, offset, corn_lat)
          call resize_nemo_field(corn_lon,corner=.true.)
          call resize_nemo_field(corn_lat,corner=.true.)

          !--- clon and clat are no longer required
          deallocate( clon, clat )

          !--- Integer mask for ocean grid cells (1=ocean, 0=land)
          if ( associated(imask) ) deallocate(imask)
          allocate( imask(nx,ny) )
          imask = nemo_tmask
          !--- Remove overlap longitudes
          call resize_nemo_field(imask)

          !--- Grid cell area
          if ( associated(area) ) deallocate(area)
          allocate( area(nx,ny) )
          area = nemo_e1t * nemo_e2t
          !--- Remove overlap longitudes
          call resize_nemo_field(area)

          !--- Since overlap longitudes were removed we must reset nx for use below
          nx=nx-2

          if ( verbose > 10 ) then
            write(6,*)"define_grid_from_nemo_coords: Boundary grid cells"
            nxc=size(corn_lon,dim=1)
            do idx1=1,size(mid_lon,dim=2)
              write(6,*)"-----------------------------------------------------------------------------"
              !--- Center at (i,j) in mid_lon and mid_lat
              !--- Corners in corn_lon and corn_lat at
              !--- SW corner at (i  ,j  )
              !--- SE corner at (i+1,j  )
              !--- NW corner at (i  ,j+1)
              !--- NE corner at (i+1,j+1)
              write(6,'("   lat index: ",i5)')idx1
              write(6,'(" cell center: left boundary ",2f12.6," right boundary ",2f12.6)') &
                  mid_lon(1,idx1),mid_lat(1,idx1),mid_lon(nxc-1,idx1),mid_lat(nxc-1,idx1)
              write(6,'("integer mask: left boundary ",f12.6,12x," right boundary ",f12.6)') &
                  imask(1,idx1),imask(nxc-1,idx1)
              write(6,'("   SW corner: left boundary ",2f12.6," right boundary ",2f12.6)') &
                  corn_lon(1,idx1),corn_lat(1,idx1),corn_lon(nxc-1,idx1),corn_lat(nxc-1,idx1)
              write(6,'("   SE corner: left boundary ",2f12.6," right boundary ",2f12.6)') &
                  corn_lon(2,idx1),corn_lat(2,idx1),corn_lon(nxc,idx1),corn_lat(nxc,idx1)
              write(6,'("   NW corner: left boundary ",2f12.6," right boundary ",2f12.6)') &
                  corn_lon(1,idx1+1),corn_lat(1,idx1+1),corn_lon(nxc-1,idx1+1),corn_lat(nxc-1,idx1+1)
              write(6,'("   NE corner: left boundary ",2f12.6," right boundary ",2f12.6)') &
                  corn_lon(2,idx1+1),corn_lat(2,idx1+1),corn_lon(nxc,idx1+1),corn_lat(nxc,idx1+1)
              call flush(6)
            enddo
          endif

        case ("u", "U")
          !--- The U grid requires
          !--- nemo_glamu
          !--- nemo_gphiu
          !--- nemo_glamv
          !--- nemo_gphiv
          !--- nemo_umask
          !--- nemo_e1u
          !--- nemo_e2u
          !--- offset = (/1,0/)
          gtype = 'U'
          offset = (/1,0/)
          write(6,*)"define_grid_from_nemo_coords: Grid type U not currently supported"
          call flush(6)
          call xit("define_grid_from_nemo_coords",-9)

          ! lon_name = "glamu"
          ! lat_name = "gphiu"
          ! lon_corn_name = "glamv"
          ! lat_corn_name = "gphiv"
          ! mask_name = "umask"
          ! e1_name = "e1u"
          ! e2_name = "e2u"

        case ("v", "V")
          !--- The V grid requires
          !--- nemo_glamv
          !--- nemo_gphiv
          !--- nemo_glamu
          !--- nemo_gphiu
          !--- nemo_vmask
          !--- nemo_e1v
          !--- nemo_e2v
          !--- offset = (/0,1/)
          gtype = 'V'
          offset = (/0,1/)
          write(6,*)"define_grid_from_nemo_coords: Grid type V not currently supported"
          call flush(6)
          call xit("define_grid_from_nemo_coords",-9)

          ! lon_name = "glamv"
          ! lat_name = "gphiv"
          ! lon_corn_name = "glamu"
          ! lat_corn_name = "gphiu"
          ! mask_name = "vmask"
          ! e1_name = "e1v"
          ! e2_name = "e2v"

        case ("f", "F")
          !--- The F grid requires
          !--- nemo_glamf
          !--- nemo_gphif
          !--- nemo_glamt
          !--- nemo_gphit
          !--- nemo_fmask
          !--- nemo_e1f
          !--- nemo_e2f
          !--- offset = (/-1,0/)

          !--- The values for lon_corn_name, lat_corn_name and the values for the offset
          !--- need to be verified for the F grid
          write(6,*)"define_grid_from_nemo_coords: Grid type F not currently supported"
          call flush(6)
          call xit("define_grid_from_nemo_coords",-9)
          gtype = 'F'
          offset = (/-1,0/)

          ! lon_name = "glamf"
          ! lat_name = "gphif"
          ! lon_corn_name = "glamt"
          ! lat_corn_name = "gphit"
          ! mask_name = "fmask"
          ! e1_name = "e1f"
          ! e2_name = "e2f"

        case default
          write(6,*)"define_grid_from_nemo_coords: Invalid grid type ",gtype
          call flush(6)
          call xit("define_grid_from_nemo_coords",-19)

      end select

      !--- At this point the following fields must be defined
      if ( .not. associated(mid_lon) ) then
        write(6,*)"define_grid_from_nemo_coords: mid_lon must be defined"
        call xit("define_grid_from_nemo_coords",-20)
      endif
      if ( .not. associated(mid_lat) ) then
        write(6,*)"define_grid_from_nemo_coords: mid_lat must be defined"
        call xit("define_grid_from_nemo_coords",-21)
      endif
      if ( .not. associated(corn_lon) ) then
        write(6,*)"define_grid_from_nemo_coords: corn_lon must be defined"
        call xit("define_grid_from_nemo_coords",-22)
      endif
      if ( .not. associated(corn_lat) ) then
        write(6,*)"define_grid_from_nemo_coords: corn_lat must be defined"
        call xit("define_grid_from_nemo_coords",-23)
      endif
      if ( .not. associated(imask) ) then
        write(6,*)"define_grid_from_nemo_coords: imask must be defined"
        call xit("define_grid_from_nemo_coords",-24)
      endif
      if ( .not. associated(area) ) then
        write(6,*)"define_grid_from_nemo_coords: area must be defined"
        call xit("define_grid_from_nemo_coords",-25)
      endif

      !--- Define regular decomposition counts for each dimension
      regDecomp(1) = cpl_PetCount
      regDecomp(2) = 1

      !--- Initialize the grid_id cpl_grid element
      fidmax = 0
      if ( n_grid_list > 0 ) then
        !--- Determine the maximum file ID assigned to this grid description file
        do idx1=1,n_grid_list
          !--- Ignore grids that do not have file_name defined
          if ( len_trim(cpl_grid(idx1)%file_name) < 1 ) cycle
          !--- Extract the file name and file ID from the current cpl_grid file name
          strng = grid_desc_file_name(idx1, fid=fid)
          !--- If the cpl_grid file name matches the current file name then assign fidmax
          if ( trim(adjustl(strng)) .eq. "internal" ) fidmax = max(fid,fidmax)
        enddo
      endif
      !--- The file ID for the current grid is one greater than the maximum file ID
      !--- that was found in cpl_grid for this file name or 1 if no grids matched
      fid = fidmax + 1
      sfx = " "
      write(sfx,'("(",i3.3,")")')fid

      cpl_grid(grid_id)%file_name = "internal"//trim(sfx)
      cpl_grid(grid_id)%title     = trim(grid_title)
      cpl_grid(grid_id)%nx        = nx
      cpl_grid(grid_id)%ny        = ny
      cpl_grid(grid_id)%size      = nx*ny
      cpl_grid(grid_id)%edges     = 0

      !--- Create a grid containing 1 periodic dimension and a regular distribution
      gridEdgeLWidth = (/ 0, 0 /)
      gridEdgeUWidth = (/ 0, 1 /)
      nx_esmf = nx
      ny_esmf = ny
      grid_list(grid_id) = ESMF_GridCreate1PeriDim( regDecomp=regDecomp, &
                                                    maxIndex=(/nx_esmf,ny_esmf/), &
                                                    coordSys=ESMF_COORDSYS_SPH_DEG, &
                                                    gridEdgeLWidth=gridEdgeLWidth,  &
                                                    gridEdgeUWidth=gridEdgeUWidth,  &
                                                    indexflag=ESMF_INDEX_GLOBAL,    &
                                                    name=trim(grid_title), &
                                                    rc=rc )
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot create grid from NEMO coordinates", line=lnum, method=method)

      !--- Allocate storage for the grid center coordinates
      call ESMF_GridAddCoord(grid_list(grid_id), staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot allocate space for grid center coords", line=lnum, method=method)

      !--- Allocate storage for the grid corner coordinates
      call ESMF_GridAddCoord(grid_list(grid_id), staggerloc=ESMF_STAGGERLOC_CORNER, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot allocate space for grid corner coords", line=lnum, method=method)

      !--- Allocate storage for the grid mask
      call ESMF_GridAddItem(grid_list(grid_id), itemflag=ESMF_GRIDITEM_MASK, &
                                            staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot allocate space for grid mask", line=lnum, method=method)

      !--- Retrieve a pointer to the grid mask
      !--- add assign these values using data read from the file
      call ESMF_GridGetItem(grid_list(grid_id), &
                            staggerloc=ESMF_STAGGERLOC_CENTER, &
                            itemflag=ESMF_GRIDITEM_MASK, &
                            farrayPtr=grid_mask, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot retrieve grid mask", line=lnum, method=method)

      nxc = size(grid_mask,dim=1)
      nyc = size(grid_mask,dim=2)
      if ( nxc /= size(imask,dim=1) .or. nyc /= size(imask,dim=2) ) then
        write(6,*)"define_grid_from_nemo_coords: Shape mismatch for grid imask"
        write(6,*)"    shape expected = ",shape(imask)
        write(6,*)"       shape found = ",nxc,nyc
        call xit("define_grid_from_nemo_coords",-1)
      endif
      grid_mask = imask(1:nxc,1:nyc)

      !--- Allocate and assign space for this mask in cpl_grid
      if ( associated(cpl_grid(grid_id)%mask) ) deallocate(cpl_grid(grid_id)%mask)
      allocate( cpl_grid(grid_id)%mask(nxc,nyc) )
      cpl_grid(grid_id)%mask = imask(1:nxc,1:nyc)

      !--- This is a spherical grid containing 1 periodic dimension (lon)
      !--- 1st coord (lon) is 2D and depends on both Grid dimensions
      !--- 2nd coord (lat) is 2D and depends on both Grid dimensions

      !--- Retrieve a pointer to the cell center lon (dim 1) coordinate values
      !--- add assign these values using coordinate data defined above
      call ESMF_GridGetCoord(grid_list(grid_id), coordDim=one, localDE=zero, &
                             staggerloc=ESMF_STAGGERLOC_CENTER, &
                             farrayPtr=grid_mid_lon, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to retrieve grid_mid_lon coord values", line=lnum, method=method)

      nxc = size(grid_mid_lon,dim=1)
      nyc = size(grid_mid_lon,dim=2)
      if ( nxc /= size(mid_lon,dim=1) .or. nyc /= size(mid_lon,dim=2) ) then
        write(6,*)"define_grid_from_nemo_coords: Shape mismatch for center longitude"
        write(6,*)"    shape expected = ",shape(mid_lon)
        write(6,*)"       shape found = ",nxc,nyc
        call xit("define_grid_from_nemo_coords",-1)
      endif
      grid_mid_lon = mid_lon(1:nxc,1:nyc)

      !--- Retrieve a pointer to the cell center lat (dim 2) coordinate values
      !--- add assign these values using coordinate data defined above
      call ESMF_GridGetCoord(grid_list(grid_id), coordDim=two, localDE=zero, &
                             staggerloc=ESMF_STAGGERLOC_CENTER, &
                             farrayPtr=grid_mid_lat, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to retrieve grid_mid_lat coord values", line=lnum, method=method)

      nxc = size(grid_mid_lat,dim=1)
      nyc = size(grid_mid_lat,dim=2)
      if ( nxc /= size(mid_lat,dim=1) .or. nyc /= size(mid_lat,dim=2) ) then
        write(6,*)"define_grid_from_nemo_coords: Shape mismatch for center latitude"
        write(6,*)"    shape expected = ",shape(mid_lat)
        write(6,*)"       shape found = ",nxc,nyc
        call xit("define_grid_from_nemo_coords",-1)
      endif
      grid_mid_lat = mid_lat(1:nxc,1:nyc)

      !--- Retrieve a pointer to the cell corner lon (dim 1) coordinate values
      !--- add assign these values using coordinate data defined above
      call ESMF_GridGetCoord(grid_list(grid_id), coordDim=one, localDE=zero, &
                             staggerloc=ESMF_STAGGERLOC_CORNER, &
                             farrayPtr=grid_corn_lon, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to retrieve grid_corn_lon coord values", line=lnum, method=method)

      !--- This grid is periodic in the first dimension (lon) so there is no need
      !--- for an extra satgger location in longitude. In fact this is not allowed
      !--- by ESMF_GridCreate1PeriDim. The second dimension (lat) does need 1 extra
      !--- stagger location. Therefore the shape of the corner stagger location
      !--- arrays should be (nx,ny+1) where (nx,ny) is the shape of the grid center
      !--- coordinates. The corn_lon and corn_lat arrays are dimensioned (nx+1,ny+1)
      !--- Thus the the check on shape of the corner stagger arrays is as follows.
      nxc = size(grid_corn_lon,dim=1)
      nyc = size(grid_corn_lon,dim=2)
      if ( nxc /= size(corn_lon,dim=1)-1 .or. nyc /= size(corn_lon,dim=2) ) then
        write(6,*)"define_grid_from_nemo_coords: Shape mismatch for corner longitude"
        write(6,*)"    shape expected = ",size(corn_lon,dim=1)-1,size(corn_lon,dim=2)
        write(6,*)"       shape found = ",nxc,nyc
        call xit("define_grid_from_nemo_coords",-1)
      endif
      grid_corn_lon = corn_lon(1:nxc,1:nyc)

      !--- Retrieve a pointer to the cell corner lat (dim 2) coordinate values
      !--- add assign these values using coordinate data defined above
      call ESMF_GridGetCoord(grid_list(grid_id), coordDim=two, localDE=zero, &
                             staggerloc=ESMF_STAGGERLOC_CORNER, &
                             farrayPtr=grid_corn_lat, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to retrieve grid_corn_lat coord values", line=lnum, method=method)

      nxc = size(grid_corn_lat,dim=1)
      nyc = size(grid_corn_lat,dim=2)
      if ( nxc /= size(corn_lat,dim=1)-1 .or. nyc /= size(corn_lat,dim=2) ) then
        write(6,*)"define_grid_from_nemo_coords: Shape mismatch for corner latitude"
        write(6,*)"    shape expected = ",size(corn_lat,dim=1)-1,size(corn_lat,dim=2)
        write(6,*)"       shape found = ",nxc,nyc
        call xit("define_grid_from_nemo_coords",-1)
      endif
      grid_corn_lat = corn_lat(1:nxc,1:nyc)

      !--- Grid cell area

      !--- Assign this cell area data to the global array ocn_grid_cell_area
      !--- Since nemo_et1 and nemo_e2t have units of meters this area will be in meters squared
      if ( associated(ocn_grid_cell_area) ) deallocate(ocn_grid_cell_area)
      allocate( ocn_grid_cell_area(size(area,dim=1), size(area,dim=2)) )
      ocn_grid_cell_area = area
      garea = sum(area)
      write(6,*)"define_grid_from_nemo_coords: Dimensional NEMO total grid area ",garea

      use_ocn_grid_cell_area = .true.
      if ( use_ocn_grid_cell_area ) then
        !--- Use the area found in the ocn_grid_cell_area array
        !--- Renormalize the area to the unit sphere for use in regridding
        if ( atm_earth_radius < 0.0 ) then
          !--- It is an error if this is not defined
          write(6,*)"define_grid_from_nemo_coords: atm_earth_radius is not defined"
          call xit("define_grid_from_nemo_coords",-1)
        endif
        !--- Use the same scaling for the ocean as for the AGCM
        rad2 = atm_earth_radius * atm_earth_radius
        !INFO: using NEMO earth radius of 6371229 meters rather than AGCM earth radius 6371220
        !INFO: rad2 = 6371229.0_8 * 6371229.0_8
        !INFO: This change results in a slight difference in the normalized sum of the total ocean area
        !INFO: from 8.94972230908983235 when earth radius is 6371229
        !INFO: to   8.94969702434640624 when earth radius is 6371220
        !INFO: For comparison, the normalized sum of the total ocean area on the AGCM grid
        !INFO: is   8.94950863406954689 when earth radius is 6371220
        !INFO: Allowing ESMF to calculate grid cell areas results in the sum of the total ocean areas
        !INFO:    8.93259740050093498 for the AGCM grid
        !INFO:    8.93259740050178053 for the NEMO grid
        area = area/rad2
        garea = sum(area)
        write(6,*)"define_grid_from_nemo_coords:  Normalized total ocean area on NEMO grid ",garea
      else
        !--- Calculate the grid cell area internally
        !--- Define a real kind=8 field using the current grid
        call compute_grid_area(grid_id, area_data1)

        if ( size(area,dim=1) /= size(area_data1,dim=1) .or. &
             size(area,dim=2) /= size(area_data1,dim=2) ) then
          write(6,*)"define_grid_from_nemo_coords: Shape mismatch for grid areas"
          write(6,*)"expected nx,ny = ",size(area,dim=1),size(area,dim=2)
          write(6,*)"   found nx,ny = ",size(area_data1,dim=1),size(area_data1,dim=2)
          call xit("define_grid_from_nemo_coords",-1)
        endif

        area = area_data1
      endif

      !--- Allocate space for the area grid item in this grid
      call ESMF_GridAddItem(grid_list(grid_id), itemflag=ESMF_GRIDITEM_AREA, &
           staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot allocate space for grid area", line=lnum, method=method)

      !--- Retrieve a pointer to the grid area in the current grid
      call ESMF_GridGetItem(grid_list(grid_id), &
                            staggerloc=ESMF_STAGGERLOC_CENTER, &
                            itemflag=ESMF_GRIDITEM_AREA, &
                            farrayPtr=area_data2, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Cannot retrieve grid area", line=lnum, method=method)

      !--- Verify that the referenced array sizes are equal
      if ( size(area,dim=1) /= size(area_data2,dim=1) .or. &
           size(area,dim=2) /= size(area_data2,dim=2) ) then
        write(6,*)"define_grid_from_nemo_coords: Shape mismatch for grid areas"
        write(6,*)"expected nx,ny = ",size(area,dim=1),size(area,dim=2)
        write(6,*)"   found nx,ny = ",size(area_data2,dim=1),size(area_data2,dim=2)
        call xit("define_grid_from_nemo_coords",-1)
      endif

      !--- Assign the grid area in the current grid with the data from area_field
      area_data2 = area

      if (verbose > -10) then
        write(6,*)"define_grid_from_nemo_coords:    Internal NEMO grid total area is ",sum(area_data2)
        call flush(6)
      endif

      !--- Clean up area arrays
      nullify(area_data1)
      nullify(area_data2)

      !--- Associate this grid with the cpl_grid array
      cpl_grid(grid_id)%grid => grid_list(grid_id)

      if ( associated(mid_lon) )  deallocate(mid_lon)
      if ( associated(mid_lat) )  deallocate(mid_lat)
      if ( associated(corn_lon) ) deallocate(corn_lon)
      if ( associated(corn_lat) ) deallocate(corn_lat)
      if ( associated(imask) )    deallocate(imask)
      if ( associated(area) )     deallocate(area)
      if ( associated(awrk2d) )   deallocate(awrk2d)
      if ( associated(bwrk2d) )   deallocate(bwrk2d)
      if ( associated(iwrk2d) )   deallocate(iwrk2d)

    end subroutine define_grid_from_nemo_coords

    !***************************************************************************
    !--- Print info about a particular grid
    !***************************************************************************
    subroutine print_grid_info(this_grid)
      type(ESMF_Grid), intent(in) :: this_grid

      !--- Local
      character(64) :: method = "print_grid_info"
      integer(ESMF_KIND_I) :: lnum

      type(ESMF_TypeKind_Flag)       :: grid_coord_kind
      type(ESMF_CoordSys_Flag)       :: coord_sys
      type(ESMF_Index_Flag)          :: grid_indexflag
      type(ESMF_GridStatus_Flag)     :: grid_status
      type(ESMF_DistGrid)            :: distgrid

      integer(ESMF_KIND_I) :: dim_count, tile_count, staggerloc_count, localde_count
      integer(ESMF_KIND_I), target :: coord_dim_count(2)
      integer(ESMF_KIND_I), target :: coord_dim_map(2,2)
      integer(ESMF_KIND_I) :: rc, idx, nx, ny, rank
      character(128) :: grid_name
      real(ESMF_KIND_R8), pointer :: coord2d(:,:) => null()
      real(ESMF_KIND_R8), pointer :: coord1d(:) => null()
      type(ESMF_Array) :: coordAry

      call ESMF_GridGet(this_grid, rc = rc, &
                          coordTypeKind = grid_coord_kind,    &
                               dimCount = dim_count,          &
                              tileCount = tile_count,         &
                        staggerlocCount = staggerloc_count,   &
                           localDECount = localde_count,      &
                               distgrid = distgrid,           &
                               coordSys = coord_sys,          &
                          coordDimCount = coord_dim_count,    &
                            coordDimMap = coord_dim_map,      &
                              indexflag = grid_indexflag,     &
                                 status = grid_status,        &
                                   name = grid_name)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Problem getting grid info", line=lnum, method=method)

      write(*,*)'print_grid_info: ESMF grid name: ',trim(grid_name)
      write(*,*)'print_grid_info:             coord system = ',coord_sys
      write(*,*)'print_grid_info:                dim count = ',dim_count
      write(*,*)'print_grid_info:          coord dim count = ',coord_dim_count
      write(*,*)'print_grid_info:          coord dim map   = ',coord_dim_map
      write(*,*)'print_grid_info:          coord type kind = ',grid_coord_kind
      write(*,*)'print_grid_info:   stagger location count = ',staggerloc_count
      write(*,*)'print_grid_info:               tile count = ',tile_count
      write(*,*)'print_grid_info:           local DE count = ',localde_count
      write(*,*)'print_grid_info:               index flag = ',grid_indexflag
      write(*,*)'print_grid_info:              status flag = ',grid_status

      if ( verbose > 10 ) then
        call ESMF_DistGridPrint(distgrid, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Problem printing DisGrid", line=lnum, method=method)
      endif

      !--- Print dimension 1 center values
      call ESMF_GridGetCoord(this_grid, coordDim=one, array=coordAry, &
                             staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to retrieve dim 1 center coord values", line=lnum, method=method)

      call ESMF_ArrayGet(coordAry, rank=rank, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to retrieve rank of dim 1 center coord values", line=lnum, method=method)
      if ( rank == 1 ) then
        call ESMF_ArrayGet(coordAry, farrayPtr=coord1d, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve 1d pointer to dim 1 center coord values", line=lnum, method=method)
        write(6,*)"print_grid_info: Dim 1 center coords = ",coord1d
        nullify(coord1d)
      else if ( rank == 2 ) then
        call ESMF_ArrayGet(coordAry, farrayPtr=coord2d, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve 2d pointer to dim 1 center coord values", line=lnum, method=method)
        write(6,*)"print_grid_info: Dim 1 center coords row = ",coord2d(1:5,1)," ..."
        write(6,*)"print_grid_info: Dim 1 center coords col = ",coord2d(1,1:5)," ..."
        nullify(coord2d)
      else
        write(*,*)'print_grid_info: Invalid rank = ',rank
        call xit("print_grid_info",-1)
      endif

      !--- Print dimension 1 corner values
      call ESMF_GridGetCoord(this_grid, coordDim=one, array=coordAry, &
                             staggerloc=ESMF_STAGGERLOC_CORNER, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to retrieve dim 1 corner coord values", line=lnum, method=method)

      call ESMF_ArrayGet(coordAry, rank=rank, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to retrieve rank of dim 1 corner coord values", line=lnum, method=method)
      if ( rank == 1 ) then
        call ESMF_ArrayGet(coordAry, farrayPtr=coord1d, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve 1d pointer to dim 1 corner coord values", line=lnum, method=method)
        write(6,*)"print_grid_info: Dim 1 corner coords = ",coord1d
        nullify(coord1d)
      else if ( rank == 2 ) then
        call ESMF_ArrayGet(coordAry, farrayPtr=coord2d, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve 2d pointer to dim 1 corner coord values", line=lnum, method=method)
        write(6,*)"print_grid_info: Dim 1 corner coords row = ",coord2d(1:5,1)," ..."
        write(6,*)"print_grid_info: Dim 1 corner coords col = ",coord2d(1,1:5)," ..."
        nullify(coord2d)
      else
        write(*,*)'print_grid_info: Invalid rank = ',rank
        call xit("print_grid_info",-1)
      endif

      !--- Print dimension 2 center values
      call ESMF_GridGetCoord(this_grid, coordDim=two, array=coordAry, &
                             staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to retrieve dim 2 center coord values", line=lnum, method=method)

      call ESMF_ArrayGet(coordAry, rank=rank, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to retrieve rank of dim 2 center coord values", line=lnum, method=method)
      if ( rank == 1 ) then
        call ESMF_ArrayGet(coordAry, farrayPtr=coord1d, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve 1d pointer to dim 2 center coord values", line=lnum, method=method)
        write(6,*)"print_grid_info: Dim 2 center coords = ",coord1d
        nullify(coord1d)
      else if ( rank == 2 ) then
        call ESMF_ArrayGet(coordAry, farrayPtr=coord2d, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve 2d pointer to dim 2 center coord values", line=lnum, method=method)
        write(6,*)"print_grid_info: Dim 2 center coords row = ",coord2d(1:5,1)," ..."
        write(6,*)"print_grid_info: Dim 2 center coords col = ",coord2d(1,1:5)," ..."
        nullify(coord2d)
      else
        write(*,*)'print_grid_info: Invalid rank = ',rank
        call xit("print_grid_info",-1)
      endif

      !--- Print dimension 2 corner values
      call ESMF_GridGetCoord(this_grid, coordDim=two, array=coordAry, &
                             staggerloc=ESMF_STAGGERLOC_CORNER, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to retrieve dim 2 corner coord values", line=lnum, method=method)

      call ESMF_ArrayGet(coordAry, rank=rank, rc=rc)
      lnum = __LINE__ ; call handle_ESMF_err(rc, &
          "Unable to retrieve rank of dim 2 corner coord values", line=lnum, method=method)
      if ( rank == 1 ) then
        call ESMF_ArrayGet(coordAry, farrayPtr=coord1d, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve 1d pointer to dim 2 corner coord values", line=lnum, method=method)
        write(6,*)"print_grid_info: Dim 2 corner coords = ",coord1d
        nullify(coord1d)
      else if ( rank == 2 ) then
        call ESMF_ArrayGet(coordAry, farrayPtr=coord2d, rc=rc)
        lnum = __LINE__ ; call handle_ESMF_err(rc, &
            "Unable to retrieve 2d pointer to dim 2 corner coord values", line=lnum, method=method)
        write(6,*)"print_grid_info: Dim 2 corner coords row = ",coord2d(1:5,1)," ..."
        write(6,*)"print_grid_info: Dim 2 corner coords col = ",coord2d(1,1:5)," ..."
        nullify(coord2d)
      else
        write(*,*)'print_grid_info: Invalid rank = ',rank
        call xit("print_grid_info",-1)
      endif

      !--- NOTE: Do not destroy the ESMF array coordAry here (do not call ESMF_ArrayDestroy)
      !---       since this call appears to deallocate coordinate values

    end subroutine print_grid_info

    !***************************************************************************
    !--- Given a storage record, return a pointer to an ESMF array containing
    !--- data from the storage record
    !***************************************************************************
    function store2Array(rec_ptr) result(array_ptr)

      implicit none

      type(store_record_t), pointer :: rec_ptr
      type(ESMF_Array), pointer :: array_ptr

      !--- Local
      type(cpl_vinfo_t) :: curr_vinfo
      integer(ESMF_KIND_I) :: rc, lnum
      character(64) :: method = "store2Array"
      type(ESMF_DistGrid)  :: distgrid
      type(ESMF_ArraySpec) :: arrayspec
      integer(ESMF_KIND_I) :: rank, regDecomp(2), min_index(2), max_index(2)

      rec_ptr => null()

      !--- Access information about this field
      curr_vinfo = find_cpl_vinfo( name=trim(adjustl(rec_ptr%name)) )

      !--- Define a regular decomposition for all 2D variables
      regDecomp(1) = 1
      regDecomp(2) = 1

      rank=2
      !--- Dimension the array to include overlap longitudes, if any
      min_index(1) = 1
      min_index(2) = 1
      max_index(1) = curr_vinfo%nlon + curr_vinfo%olap
      max_index(2) = curr_vinfo%nlat

      call ESMF_ArraySpecSet(arrayspec, typekind=ESMF_TYPEKIND_R8, rank=rank, rc=rc)

      distgrid = ESMF_DistGridCreate( min_index, max_index, regDecomp=regDecomp, rc=rc)

      allocate( array_ptr )
      array_ptr = ESMF_ArrayCreate(distgrid, arrayspec, name=trim(rec_ptr%name), rc=rc)

      call ESMF_DistGridDestroy(distgrid, rc=rc)

    end function store2Array

end module ccc_esmf
