module grid_desc

  use netcdf
  use ccc_io

  implicit none

  private

  public :: create_grid_desc
  public :: grid_desc_list, n_grid_desc_list
  public :: read_xxx_grid_desc_file
  public :: grid_desc_t

  integer, parameter ::                &
      ws_i4  = selected_int_kind(6)   ,&
      ws_i8  = selected_int_kind(13)  ,&
      ws_r4  = selected_real_kind(6)  ,&
      ws_r8  = selected_real_kind(13) ,&
      ws_r16 = selected_real_kind(26)

  real (kind=ws_r8), parameter ::  &
      zero   = 0.0_ws_r8, &
      one    = 1.0_ws_r8, &
      two    = 2.0_ws_r8, &
      three  = 3.0_ws_r8, &
      four   = 4.0_ws_r8, &
      five   = 5.0_ws_r8, &
      half   = 0.5_ws_r8, &
      quart  = 0.25_ws_r8, &
      bignum = 1.e+20_ws_r8, &
      tiny   = 1.e-14_ws_r8, &
      pi     = 3.141592653589793_ws_r8, &
      pi2    = two*pi, &
      pih    = half*pi

  integer (kind=ws_i4), parameter :: one_i4=1

  !--- A grid description is normally stored in a netcdf file containing information
  !--- about the grid cells and the size and shape of the grid.
  !--- This file may be read and stored in grid_desc_t containers using read_xxx_grid_desc_file
  !--- The grid description file may also be created using the subroutine create_grid_desc
  !--- The grid description file must be in a particular format similar to the following
  !---    dimensions:
  !---            grid_size = 8192 ;
  !---            grid_corners = 4 ;
  !---            grid_rank = 2 ;
  !---    variables:
  !---            int grid_dims(grid_rank) ;
  !---            double grid_center_lat(grid_size) ;
  !---                    grid_center_lat:units = "degrees" ;
  !---            double grid_center_lon(grid_size) ;
  !---                    grid_center_lon:units = "degrees" ;
  !---            int grid_imask(grid_size) ;
  !---                    grid_imask:units = "unitless" ;
  !--- NOTE: Netcdf files contain multi-dimensional arrays are ordered with the last dimension
  !---       varying the fastest. Associated fortran arrays should be dimensioned in reverse
  !---       order to the order in the netcdf files.
  !---       e.g. fortran grid_corner_lat|lon are dimensioned (grid_corners, grid_size)
  !---            double grid_corner_lat(grid_size, grid_corners) ;
  !---                    grid_corner_lat:units = "degrees" ;
  !---            double grid_corner_lon(grid_size, grid_corners) ;
  !---                    grid_corner_lon:units = "degrees" ;
  !---            double grid_area(grid_size) ;
  !---                    grid_area:units = "unitless" ;
  !---
  !---    // global attributes:
  !---                    :title = "T63 Gaussian Grid -- 128x64" ;

  !--- Define a container to hold a grid description
  type grid_desc_t
    !--- A name that may be used to identify this grid
    character(len=128) :: title = " "

    !--- The name of the disk file from which this file was read, if any
    character(len=128) :: file = " "

    !--- Total number of grid cells in this grid
    integer (ws_i4) :: size    = -1

    !--- Number of corners in each grid cell
    integer (ws_i4) :: corners = -1

    !--- Number of dimensions in this grid
    integer (ws_i4) :: rank    = -1

    !--- Grid coordinate units (either degrees or radians)
    character(len=16)  :: units = "degrees         "

    !--- The shape of this grid ...dimensioned (rank)
    !--- e.g. for a 2D lat/lon grid grid_dims(1:2) = (nlon, nlat)
    integer (ws_i4), dimension(:), pointer :: dims => null()

    !--- Integer mask (0 or 1) used to flag which cells participate in the mapping
    !--- dimensioned (size)
    integer (ws_i4), dimension(:), pointer :: imask => null()

    !--- Grid center latitude and longitude ...dimensioned (size)
    real (ws_r8), dimension(:), pointer :: center_lat => null()
    real (ws_r8), dimension(:), pointer :: center_lon => null()

    !--- Grid corner latitude and longitude ...dimensioned (corners, size)
    real (ws_r8), dimension(:,:), pointer :: corner_lat => null()
    real (ws_r8), dimension(:,:), pointer :: corner_lon => null()

    !--- Area of each grid cell ...dimensioned (size)
    real (ws_r8), dimension(:), pointer :: area => null()
  end type grid_desc_t

  !--- The maximum number of grid descriptions that may be held in memory
  integer, parameter :: max_grid_desc=20

  !--- The number of grid descriptions currently in use
  integer, save :: n_grid_desc_list = 0

  !--- An array containing a list of all grid descriptions currently in use
  type(grid_desc_t), save, target, dimension(max_grid_desc) :: grid_desc_list

  type grid_array_t
    type(grid_desc_t), pointer :: grid => null()
  end type grid_array_t

  !--- Variables common to many subroutines in this module

  !--- Dimension sizes written to the netcdf file containing the grid description
  integer (ws_i4) :: nx, ny, grid_size, grid_rank, grid_corners

  !--- A descriptive name for the grid
  character(128) :: grid_name =   " "

  !--- The name of the disk file from which the grid data was read, if any
  character(128) ::  grid_disk_file =   " "

  !--- The shape of the grid (e.g. (nx,ny) )
  integer (kind=ws_i4), allocatable :: grid_dims(:)

  !--- The number of overlapping grid points in each dimension
  integer (kind=ws_i4), allocatable :: grid_edges(:,:)

  !--- Integer land mask
  !---   0 means do not include this point
  !---   1 means include this point
  integer (ws_i4), allocatable  :: grid_imask(:)

  !--- lat/lon coordinates for each grid center in degrees
  real (ws_r8), allocatable :: grid_center_lat(:), grid_center_lon(:)

  !--- lat/lon coordinates for each grid cell corner in degrees
  real (ws_r8), allocatable :: grid_corner_lat(:,:), grid_corner_lon(:,:)

  !--- An optional array containing grid cell areas
  !--- Scrip routines will determine these values if not supplied by the user
  real (ws_r8), allocatable :: grid_area(:)

  !--- nav_lon, nav_lat are only required for an ORCA grid
  real(ws_r8), allocatable, dimension(:,:) :: nav_lon, nav_lat

  !--- Logical flags to control reading and/or writing grid area to files
  logical ::  read_grid_area = .false.
  logical :: write_grid_area = .false.

  !--- An integer flag used to determine the type of points to be masked
  !---   curr_mask_type = -1 ..no mask (all points particpate in remapping)
  !---   curr_mask_type =  0 ..use a default mask depending on the type of grid
  !---   curr_mask_type =  1 ..mask land points  (land points are not used)
  !---   curr_mask_type =  2 ..mask ocean points (ocean points are not used)
  integer, save :: curr_mask_type = 0

  !--- usr_mask_gg is the name of the mask file used by create_grid_file_gg
  !--- If defined (non-blank) then this is the name of a CCCma binary file containing
  !--- a single record (named either MASK or GC) defining the gaussian grid land mask.
  !--- Similar (or the same) files are read by the AGCM at run time and used to create
  !--- the MASK record that is found in the "an" file. In this case the file name is
  !--- supplied by the parmsub parameter "maskfile", is processed by initg12 in the
  !--- deck inphys16.cdk (or the current version thereof) and written to the "an" file.
  !--- In the file usr_mask_gg, land = -1 and all other points are considerer not land.
!--- TODO ---
  character (128) :: usr_mask_gg = " "

  !--- usr_ocn_desc is the name of a file containing the ocean data description
  !--- This file must include the records: LON, LONH, LAT, LATH, Z, ZH, ZBOT
  !--- It may contain other records, which will be ignored
!--- TODO ---
  character (128) :: usr_ocn_desc = " "

  !-----------------------------------------------------------------------
  !--- CCCma specific variables
  !-----------------------------------------------------------------------

  !--- CCCma I/O buffer
  integer, parameter :: maxx=2099200  !--- 2048x1025
  integer :: ibuf(8),idat(2*maxx)
  common /icom/ ibuf,idat

  real :: wrk(maxx)
  logical :: ok
  logical :: remote
  integer :: nc4to8

  real(kind=ws_r8), parameter :: spval = 1e38_ws_r8
  private :: spval

contains

  !***********************************************************************
  !--- Assign element of a variable of type grid_desc_t using variable values
  !--- defined above in certain grid_* variables
  !***********************************************************************
  subroutine assign_one_grid(grid)
    type(grid_desc_t) :: grid
    grid % title   = " "
    grid % title   = trim(grid_name)
    grid % file   = " "
    if ( len_trim(grid_disk_file)>0 ) then
      grid % file   = trim(grid_disk_file)
    endif
    grid % size    = grid_size
    grid % corners = grid_corners
    grid % rank    = grid_rank
    grid % units   = " "
    grid % units   = "degrees"

    !--- Assign grid % dims
    if ( associated( grid%dims ) ) nullify( grid%dims )
    if ( allocated(grid_dims) ) then
      if ( size(grid_dims) /= grid_rank ) then
        write(6,*)'assign_one_grid: grid_dims must be of length ',grid_rank
        call xit("assign_one_grid",-1)
      endif
      allocate( grid%dims(grid_rank) )
      grid % dims(1:grid_rank) = grid_dims(1:grid_rank)
    else
      write(6,*)'assign_one_grid: grid % dims was not assigned.'
    endif

    !--- Assign grid % imask
    if ( associated( grid%imask ) ) nullify( grid%imask )
    if ( allocated(grid_imask) ) then
      if ( size(grid_imask) /= grid_size ) then
        write(6,*)'assign_one_grid: grid_imask must be of length ',grid_size
        call xit("assign_one_grid",-2)
      endif
      allocate( grid%imask(grid_size) )
      grid % imask(1:grid_size) = grid_imask(1:grid_size)
    else
      write(6,*)'assign_one_grid: grid % imask was not assigned.'
    endif

    !--- Assign grid % center_lat
    if ( associated( grid%center_lat ) ) nullify( grid%center_lat )
    if ( allocated(grid_center_lat) ) then
      if ( size(grid_center_lat) /= grid_size ) then
        write(6,*)'assign_one_grid: grid_center_lat must be of length ',grid_size
        call xit("assign_one_grid",-3)
      endif
      allocate( grid%center_lat(grid_size) )
      grid % center_lat(1:grid_size) = grid_center_lat(1:grid_size)
    else
      write(6,*)'assign_one_grid: grid % center_lat was not assigned.'
    endif

    !--- Assign grid % center_lon
    if ( associated( grid%center_lon ) ) nullify( grid%center_lon )
    if ( allocated(grid_center_lon) ) then
      if ( size(grid_center_lon) /= grid_size ) then
        write(6,*)'assign_one_grid: grid_center_lon must be of length ',grid_size
        call xit("assign_one_grid",-4)
      endif
      allocate( grid%center_lon(grid_size) )
      grid % center_lon(1:grid_size) = grid_center_lon(1:grid_size)
    else
      write(6,*)'assign_one_grid: grid % center_lon was not assigned.'
    endif

    !--- Assign grid % corner_lat
    if ( associated( grid%corner_lat ) ) nullify( grid%corner_lat )
    if ( allocated(grid_corner_lat) ) then
      if ( size(grid_corner_lat, dim=1) /= grid_corners .or. &
           size(grid_corner_lat, dim=2) /= grid_size ) then
        write(6,*)'assign_one_grid: grid_corner_lat must be of shape ', &
                   grid_corners,grid_size
        call xit("assign_one_grid",-5)
      endif
      allocate( grid%corner_lat(grid_corners,grid_size) )
      grid % corner_lat(1:grid_corners,1:grid_size) = &
        grid_corner_lat(1:grid_corners,1:grid_size)
    else
      write(6,*)'assign_one_grid: grid % corner_lat was not assigned.'
    endif

    !--- Assign grid % corner_lon
    if ( associated( grid%corner_lon ) ) nullify( grid%corner_lon )
    if ( allocated(grid_corner_lon) ) then
      if ( size(grid_corner_lon, dim=1) /= grid_corners .or. &
           size(grid_corner_lon, dim=2) /= grid_size ) then
        write(6,*)'assign_one_grid: grid_corner_lon must be of shape ', &
                   grid_corners,grid_size
        call xit("assign_one_grid",-5)
      endif
      allocate( grid%corner_lon(grid_corners,grid_size) )
      grid % corner_lon(1:grid_corners,1:grid_size) = &
        grid_corner_lon(1:grid_corners,1:grid_size)
    else
      write(6,*)'assign_one_grid: grid % corner_lon was not assigned.'
    endif

    !--- Assign grid % area
    if ( associated( grid%area ) ) nullify( grid%area )
    if ( allocated(grid_area) ) then
      if ( size(grid_area) /= grid_size ) then
        write(6,*)'assign_one_grid: grid_area must be of length ',grid_size
        call xit("assign_one_grid",-6)
      endif
      allocate( grid%area(grid_size) )
      grid % area(1:grid_size) = grid_area(1:grid_size)
    else
      write(6,*)'assign_one_grid: grid % area was not assigned.'
    endif
  end subroutine assign_one_grid

  !***********************************************************************
  !--- Print grid coordinates to stdout
  !***********************************************************************
  subroutine print_grid()
    implicit none
    integer(kind=ws_i4) :: ix, i, j
    character(len=32) :: coord_fmt
    character(len=256) :: afmt, bfmt

    coord_fmt = '"(",f8.3,",",f8.3,")"'
    afmt = '(i6,2i4,3(3x,' //trim(coord_fmt)// '))'
    bfmt = '(17x,' //trim(coord_fmt)// ',25x,' //trim(coord_fmt)// ')'
    j = 0
    do ix=1,grid_size
      if ( mod(ix,nx) == 1 ) j = j + 1
      i = ix - (j-1)*nx
      write(6,trim(afmt)) ix, i, j, &
              grid_corner_lon(4,ix), grid_corner_lat(4,ix), &
              grid_center_lon(ix),   grid_center_lat(ix),   &
              grid_corner_lon(3,ix), grid_corner_lat(3,ix)
      write(6,trim(bfmt)) &
              grid_corner_lon(1,ix), grid_corner_lat(1,ix), &
              grid_corner_lon(2,ix), grid_corner_lat(2,ix)
    enddo

  end subroutine print_grid

  subroutine print_one_cell(cent_lon,cent_lat,corn_lon,corn_lat,i,j)
    !--- lat/lon coordinates for each grid cell corner in degrees
    real (ws_r8) :: corn_lon(:,:,:), corn_lat(:,:,:)
    real (ws_r8) :: cent_lon(:,:), cent_lat(:,:)
    integer :: i, j, ix
    character(len=32) :: coord_fmt
    character(len=256) :: afmt, bfmt

    ix=i*j
    coord_fmt = '"(",f8.3,",",f8.3,")"'
    afmt = '(i6,2i4,3(3x,' //trim(coord_fmt)// '))'
    bfmt = '(17x,' //trim(coord_fmt)// ',25x,' //trim(coord_fmt)// ')'
    write(6,trim(afmt)) ix, i, j, &
            corn_lon(4,i,j), corn_lat(4,i,j), &
              cent_lon(i,j),   cent_lat(i,j), &
            corn_lon(3,i,j), corn_lat(3,i,j)
    write(6,trim(bfmt)) &
            corn_lon(1,i,j), corn_lat(1,i,j), &
            corn_lon(2,i,j), corn_lat(2,i,j)
  end subroutine print_one_cell

  !***********************************************************************
  !--- Create a grid description file
  !***********************************************************************
  subroutine create_grid_desc(grid_alias, grid_file_name, curr_grid, &
                          mask_type, mask_gg, ocn_desc, nemo_coords, olaplon, &
                          xmin, xmax, ymin, ymax)

    character(*), intent(in)           :: grid_alias
    character(*), intent(out)          :: grid_file_name
    integer,      intent(out)          :: curr_grid
    integer,      intent(in), optional :: mask_type
    character(*), intent(in), optional :: mask_gg
    character(*), intent(in), optional :: ocn_desc
    character(*), intent(in), optional :: nemo_coords
    logical,      intent(in), optional :: olaplon
    real (kind=8), intent(in), optional :: xmin, xmax, ymin, ymax

    !--- Local
    logical :: localolaplon
    real (kind=8) :: local_xmin, local_xmax, local_ymin, local_ymax

    if ( present(mask_type) ) then
      curr_mask_type = mask_type
    else
      curr_mask_type = 0
    endif

    if ( present(olaplon) ) then
      localolaplon = olaplon
    else
      localolaplon = .false.
    endif

    !--- The use of curr_grid is depreciated
    curr_grid = 0

    if ( grid_alias(1:4) .eq. "test" ) then
      !--- This is a test grid
      if ( present(xmin) ) then
        local_xmin = xmin
      else
        local_xmin = 0.0_8
      endif
      if ( present(xmax) ) then
        local_xmax = xmax
      else
        local_xmax = 6.0_8
      endif
      if ( present(ymin) ) then
        local_ymin = ymin
      else
        local_ymin = 0.0_8
      endif
      if ( present(ymax) ) then
        local_ymax = ymax
      else
        local_ymax = 6.0_8
      endif
      call create_grid_file_test( trim(adjustl(grid_alias)), grid_file_name, &
           xmin=local_xmin, xmax=local_xmax, ymin=local_ymin, ymax=local_ymax )

    else
      !--- Call the appropriate subroutine (determined by input parameters)
      select case ( trim(adjustl(grid_alias)) )
        case ( 't31', 'T31', 't47', 'T47', 't63', 'T63', 't95', 'T95', 't127', 'T127' )
          !--- This is a gaussian grid
          if ( present(mask_gg) ) then
            call create_grid_file_gg( trim(adjustl(grid_alias)), grid_file_name, curr_grid, &
                                      localolaplon, mask_gg=trim(mask_gg) )
          else
            call create_grid_file_gg( trim(adjustl(grid_alias)), grid_file_name, curr_grid, &
                                      localolaplon )
          endif

        case ( 'ocnt', 'ocnt_256x192', 'ocnv', 'ocnv_256x192', 'use_ocn_desc' )
          !--- This is a CCCma ocean grid
          if ( present(ocn_desc) ) then
            call create_grid_file_ocn( trim(adjustl(grid_alias)), grid_file_name, curr_grid, &
                                       localolaplon, ocn_desc=trim(ocn_desc) )
          else
            call create_grid_file_ocn( trim(adjustl(grid_alias)), grid_file_name, curr_grid, &
                                       localolaplon )
          endif

        case ( 'orca2', 'orca1', 'orca025' )
          !--- This is an ORCA ocean grid
          if ( present(nemo_coords) ) then
            call create_grid_file_orca( trim(adjustl(grid_alias)), grid_file_name, curr_grid, &
                                        localolaplon, nemo_coords=trim(nemo_coords) )
          else
            call create_grid_file_orca( trim(adjustl(grid_alias)), grid_file_name, curr_grid, &
                                        localolaplon )
          endif

        case default
          write(6,*)'create_grid_desc: Invalid grid description ',trim(grid_alias)
          call xit("create_grid_desc",-1)

      end select
    endif

  end subroutine create_grid_desc

  !***********************************************************************
  !--- Create a grid description file for a test grid
  !***********************************************************************
  subroutine create_grid_file_test(grid_alias, grid_file_out, xmin, xmax, ymin, ymax)

    character(*), intent(in)            :: grid_alias
    character(*), intent(inout)         :: grid_file_out
    real (kind=8), intent(in), optional :: xmin, xmax, ymin, ymax

    !--- Local
    real (kind=8) :: dx, dy, min_v, max_v, center_v
    real (kind=8) :: local_xmin, local_xmax, local_ymin, local_ymax
    real (kind=8) :: xoffset, yoffset
    integer(kind=8) :: ix, iy, idx1, idx2
    character(256) :: local_grid_alias
    character(256) :: strng
    integer :: offset

    local_grid_alias = trim(adjustl(grid_alias))

    !--- Define the name of the output grid file
    grid_file_out = 'grid_desc_'//trim(adjustl(grid_alias))//'.nc'

    !--- Define a default x/y range here
    !--- these values may be modified via optional parameters
    if ( present(xmin) ) then
      local_xmin = xmin
    else
      local_xmin = 0.0_8
    endif
    if ( present(xmax) ) then
      local_xmax = xmax
    else
      local_xmax = 6.0_8
    endif
    if ( present(ymin) ) then
      local_ymin = ymin
    else
      local_ymin = 0.0_8
    endif
    if ( present(ymax) ) then
      local_ymax = ymax
    else
      local_ymax = 6.0_8
    endif

    xoffset = 1.0_8
    yoffset = 1.0_8

    if ( local_grid_alias(5:10) .eq. "offset" ) then
      strng(1:4) = local_grid_alias(1:4)
      strng(5:) = trim(local_grid_alias(11:))
      local_grid_alias = trim(strng)
      offset = 1
      local_xmin = local_xmin + xoffset
      local_xmax = local_xmax + xoffset
      local_ymin = local_ymin + yoffset
      local_ymax = local_ymax + yoffset
    else if ( local_grid_alias(5:11) .eq. "xoffset" ) then
      strng(1:4) = local_grid_alias(1:4)
      strng(5:) = trim(local_grid_alias(12:))
      local_grid_alias = trim(strng)
      offset = 2
      local_xmin = local_xmin + xoffset
      local_xmax = local_xmax + xoffset
    else if ( local_grid_alias(5:11) .eq. "yoffset" ) then
      strng(1:4) = local_grid_alias(1:4)
      strng(5:) = trim(local_grid_alias(12:))
      local_grid_alias = trim(strng)
      offset = 3
      local_ymin = local_ymin + yoffset
      local_ymax = local_ymax + yoffset
    else
      offset = 0
    endif

    write(6,*)"create_grid_file_test: offset=",offset,"  grid_alias=",trim(local_grid_alias)

    select case ( trim(adjustl(local_grid_alias)) )
      case ( 'test1x1' )
        !--- Regular 1x1 cartesian grid
        nx =  1
        ny =  1
        grid_name = 'Test cartesian grid -- 1x1'
      case ( 'test2x2' )
        !--- Regular 2x2 cartesian grid
        nx =  2
        ny =  2
        grid_name = 'Test cartesian grid -- 2x2'
      case ( 'test3x3' )
        !--- Regular 3x3 cartesian grid
        nx =  3
        ny =  3
        grid_name = 'Test cartesian grid -- 3x3'
      case ( 'test4x4' )
        !--- Regular 4x4 cartesian grid
        nx =  4
        ny =  4
        grid_name = 'Test cartesian grid -- 4x4'
      case ( 'test5x5' )
        !--- Regular 5x5 cartesian grid
        nx =  5
        ny =  5
        grid_name = 'Test cartesian grid -- 5x5'
      case ( 'test6x6' )
        !--- Regular 3x3 cartesian grid
        nx =  6
        ny =  6
        grid_name = 'Test cartesian grid -- 6x6'
      case ( 'test2x1' )
        nx =  2
        ny =  1
        grid_name = 'Test cartesian grid -- 2x1'
      case ( 'test1x2' )
        nx =  1
        ny =  2
        grid_name = 'Test cartesian grid -- 1x2'
      case ( 'test3x1' )
        nx =  3
        ny =  1
        grid_name = 'Test cartesian grid -- 3x1'
      case ( 'test3x2' )
        nx =  3
        ny =  2
        grid_name = 'Test cartesian grid -- 3x2'
      case ( 'test1x3' )
        nx =  1
        ny =  3
        grid_name = 'Test cartesian grid -- 1x3'
      case ( 'test2x3' )
        nx =  2
        ny =  3
        grid_name = 'Test cartesian grid -- 2x3'
      case ( 'test4x1' )
        nx =  4
        ny =  1
        grid_name = 'Test cartesian grid -- 4x1'
      case ( 'test4x2' )
        nx =  4
        ny =  2
        grid_name = 'Test cartesian grid -- 4x2'
      case ( 'test4x3' )
        nx =  4
        ny =  3
        grid_name = 'Test cartesian grid -- 4x3'
      case ( 'test1x4' )
        nx =  1
        ny =  4
        grid_name = 'Test cartesian grid -- 1x4'
      case ( 'test2x4' )
        nx =  2
        ny =  4
        grid_name = 'Test cartesian grid -- 2x4'
      case ( 'test3x4' )
        nx =  3
        ny =  4
        grid_name = 'Test cartesian grid -- 3x4'
      case ( 'test5x1' )
        nx =  5
        ny =  1
        grid_name = 'Test cartesian grid -- 5x1'
      case ( 'test5x2' )
        nx =  5
        ny =  2
        grid_name = 'Test cartesian grid -- 5x2'
      case ( 'test5x3' )
        nx =  5
        ny =  3
        grid_name = 'Test cartesian grid -- 5x3'
      case ( 'test5x4' )
        nx =  5
        ny =  4
        grid_name = 'Test cartesian grid -- 5x4'
      case ( 'test1x5' )
        nx =  1
        ny =  5
        grid_name = 'Test cartesian grid -- 1x5'
      case ( 'test2x5' )
        nx =  2
        ny =  5
        grid_name = 'Test cartesian grid -- 2x5'
      case ( 'test3x5' )
        nx =  3
        ny =  5
        grid_name = 'Test cartesian grid -- 3x5'
      case ( 'test4x5' )
        nx =  4
        ny =  5
        grid_name = 'Test cartesian grid -- 4x5'
      case ( 'test6x1' )
        nx =  6
        ny =  1
        grid_name = 'Test cartesian grid -- 6x1'
      case ( 'test6x2' )
        nx =  6
        ny =  2
        grid_name = 'Test cartesian grid -- 6x2'
      case ( 'test6x3' )
        nx =  6
        ny =  3
        grid_name = 'Test cartesian grid -- 6x3'
      case ( 'test6x4' )
        nx =  6
        ny =  4
        grid_name = 'Test cartesian grid -- 6x4'
      case ( 'test6x5' )
        nx =  6
        ny =  5
        grid_name = 'Test cartesian grid -- 6x5'
      case ( 'test1x6' )
        nx =  1
        ny =  6
        grid_name = 'Test cartesian grid -- 1x6'
      case ( 'test2x6' )
        nx =  2
        ny =  6
        grid_name = 'Test cartesian grid -- 2x6'
      case ( 'test3x6' )
        nx =  3
        ny =  6
        grid_name = 'Test cartesian grid -- 3x6'
      case ( 'test4x6' )
        nx =  4
        ny =  6
        grid_name = 'Test cartesian grid -- 4x6'
      case ( 'test5x6' )
        nx =  5
        ny =  6
        grid_name = 'Test cartesian grid -- 5x6'
      case ( 'test100x100' )
        nx =  100
        ny =  100
        grid_name = 'Test cartesian grid -- 100x100'
      case default
        write(6,'(3a)')'create_grid_file_test: Invalid grid type --> ',trim(local_grid_alias),' <--'
        call xit("create_grid_file_test",-1)
    end select

    !--- Allocate space for local arrays
    grid_size = nx*ny
    grid_rank = 2
    grid_corners = 4

    allocate( grid_dims(grid_rank) )
    allocate( grid_edges(2,grid_rank) )
    allocate( grid_imask(grid_size) )
    allocate( grid_center_lat(grid_size) )
    allocate( grid_center_lon(grid_size) )
    allocate( grid_corner_lat(grid_corners,grid_size) )
    allocate( grid_corner_lon(grid_corners,grid_size) )
    allocate( grid_area(grid_size) )

    !--- Define the grid dims array written to the grid description file
    grid_dims(1) = nx
    grid_dims(2) = ny

    !--- grid_edges is identically zero since there are no overlap cells
    grid_edges(:,:) = 0

    !-----------------------------------------------------------------------
    !--- Define longitudes (x-axis) of cell centers and cell corners
    !-----------------------------------------------------------------------
    dx = (local_xmax-local_xmin)/dble(nx)
    do iy=1,ny
      do ix=1,nx
        idx1 = (iy-1)*nx + ix

        center_v = local_xmin + dble(ix-1)*dx + half*dx
        min_v = center_v - half*dx
        max_v = center_v + half*dx
        if ( abs(min_v) < 1e-10 ) min_v = 0.0_8
        if ( abs(max_v) < 1e-10 ) max_v = 0.0_8

        grid_center_lon(idx1  ) = center_v

        !--- x value for grid corners
        grid_corner_lon(1,idx1) = min_v
        grid_corner_lon(2,idx1) = max_v
        grid_corner_lon(3,idx1) = max_v
        grid_corner_lon(4,idx1) = min_v
      end do
    end do

    !-----------------------------------------------------------------------
    !--- Define latitudes (y-axis) of cell centers and cell corners
    !-----------------------------------------------------------------------
    dy = (local_ymax-local_ymin)/dble(ny)
    do iy=1,ny
      do ix=1,nx
        idx1 = (iy-1)*nx + ix

        center_v = local_ymin + dble(iy-1)*dy + half*dy
        min_v = center_v - half*dy
        max_v = center_v + half*dy
        if ( abs(min_v) < 1e-10 ) min_v = 0.0_8
        if ( abs(max_v) < 1e-10 ) max_v = 0.0_8

        grid_center_lat(idx1  ) = center_v

        !--- y value for grid corners
        grid_corner_lat(1,idx1) = min_v
        grid_corner_lat(2,idx1) = min_v
        grid_corner_lat(3,idx1) = max_v
        grid_corner_lat(4,idx1) = max_v
      end do
    end do

    !-----------------------------------------------------------------------
    !--- grid_imask = 1 means the grid cell will be used in the remapping
    !--- grid_imask = 0 means the grid cell will not be used in the remapping
    !-----------------------------------------------------------------------
    grid_imask(:) = 1

    !-----------------------------------------------------------------------
    !--- Print this test grid
    !-----------------------------------------------------------------------
    ! call print_grid()

    !-----------------------------------------------------------------------
    !--- Create a file containing the current grid description
    !-----------------------------------------------------------------------
    call write_grid_desc_file( grid_file_out )

    !-----------------------------------------------------------------------
    !--- Clean up
    !-----------------------------------------------------------------------
    deallocate( grid_dims )
    deallocate( grid_edges )
    deallocate( grid_imask )
    deallocate( grid_center_lat )
    deallocate( grid_center_lon )
    deallocate( grid_corner_lat )
    deallocate( grid_corner_lon )
    deallocate( grid_area )

  end subroutine create_grid_file_test

  !***********************************************************************
  !--- Create a grid description file for a gaussian grid
  !***********************************************************************
  subroutine create_grid_file_gg(grid_alias, grid_file_out, curr_grid, olaplon, mask_gg)

    character(*), intent(in)           :: grid_alias
    character(*), intent(inout)        :: grid_file_out
    integer,      intent(out)          :: curr_grid
    logical,      intent(in)           :: olaplon
    character(*), intent(in), optional :: mask_gg

    !--- Local
    character (128) :: mask_gg_name = " "

    integer (kind=ws_i4) :: atm_add

    real (kind=ws_r8) :: dlon, minlon, maxlon, centerlon, minlat, maxlat, centerlat

    integer :: ix,iy,idx1,idx2
    integer :: verbose = 0

    logical :: arithmetic_mean_gg_bnds, geometric_mean_gg_bnds

    !--- Used to compute gaussian latitudes using gaussg/trigl
    real(kind=8) :: sinl(1025), wgt(1025), cosl(1025), rad(1025), wocs(1025)
    real(kind=8) :: glat(1025), glat_bnds(1026), cfac, total_area

    !-----------------------------------------------------------------------
    ! Start executable statements
    !-----------------------------------------------------------------------

    !--- If the user has supplied the name of a mask file then define it here
    !---
    !--- Note: The land mask files are CCCma binary files containing a single record
    !--- (usually named either MASK or GC). These are files that are read by the model
    !--- at run time and used to create the MASK record that is found in the "an" file.
    !--- The file name is determined at run time by the parmsub parameter "maskfile"
    !--- This file is processed by initg12 (or the current version thereof) and written
    !--- to the "an" file as part of an interm file named "ground" create in inphys16.cdk
    if ( present(mask_gg) ) then
      if (len_trim(mask_gg)>0 ) then
        mask_gg_name = trim(mask_gg)
      else
        write(6,*)'create_grid_file_gg: The land mask file name is blank.'
        call xit("create_grid_file_gg",-1)
      endif
    endif

    grid_name = " "

    select case ( trim(adjustl(grid_alias)) )
      case ( 't31', 'T31' ) !--- T31
        if (olaplon) then
          nx =  65
          grid_name = 'T31 Gaussian Grid -- 65x32'
        else
          nx =  64
          grid_name = 'T31 Gaussian Grid -- 64x32'
        endif
        ny =  32
        if ( .not. present(mask_gg) ) then
          !--- The user has not supplied the name of a file containing
          !--- the land mask. Define one here
          mask_gg_name = "gcm3.5_landmask_64x32v4"
        endif
        !--- Set the output file name here
        if (olaplon) then
          if ( curr_mask_type == -1 ) then
            grid_file_out = 'grid_desc_t31_o_nomask.nc'
          else
            grid_file_out = 'grid_desc_t31_o_'//trim(mask_gg_name)//'.nc'
          endif
        else
          if ( curr_mask_type == -1 ) then
            grid_file_out = 'grid_desc_t31_nomask.nc'
          else
            grid_file_out = 'grid_desc_t31_'//trim(mask_gg_name)//'.nc'
          endif
        endif

      case ( 't47', 'T47' ) !--- T47
        if (olaplon) then
          nx =  97
          grid_name = 'T47 Gaussian Grid -- 97x48'
        else
          nx =  96
          grid_name = 'T47 Gaussian Grid -- 96x48'
        endif
        ny =  48
        if ( .not. present(mask_gg) ) then
          !--- The user has not supplied the name of a file containing
          !--- the land mask. Define one here
          mask_gg_name = "gcm3_landmask_96x48v5"
        endif
        !--- Set the output file name here
        if (olaplon) then
          if ( curr_mask_type == -1 ) then
            grid_file_out = 'grid_desc_t47_o_nomask.nc'
          else
            grid_file_out = 'grid_desc_t47_o_'//trim(mask_gg_name)//'.nc'
          endif
        else
          if ( curr_mask_type == -1 ) then
            grid_file_out = 'grid_desc_t47_nomask.nc'
          else
            grid_file_out = 'grid_desc_t47_'//trim(mask_gg_name)//'.nc'
          endif
        endif

      case ( 't63', 'T63' ) !--- T63
        if (olaplon) then
          nx =  129
          grid_name = 'T63 Gaussian Grid -- 129x64'
        else
          nx =  128
          grid_name = 'T63 Gaussian Grid -- 128x64'
        endif
        ny =  64
        if ( .not. present(mask_gg) ) then
          !--- The user has not supplied the name of a file containing
          !--- the land mask. Define one here
          mask_gg_name = "gcm3.5_landmask_128x64v5"
        endif
        !--- Set the output file name here
        if (olaplon) then
          if ( curr_mask_type == -1 ) then
            grid_file_out = 'grid_desc_t63_o_nomask.nc'
          else
            grid_file_out = 'grid_desc_t63_o_'//trim(mask_gg_name)//'.nc'
          endif
        else
          if ( curr_mask_type == -1 ) then
            grid_file_out = 'grid_desc_t63_nomask.nc'
          else
            grid_file_out = 'grid_desc_t63_'//trim(mask_gg_name)//'.nc'
          endif
        endif

      case ( 't95', 'T95' ) !--- T95
        if (olaplon) then
          nx =  193
          grid_name = 'T95 Gaussian Grid -- 193x96'
        else
          nx =  192
          grid_name = 'T95 Gaussian Grid -- 192x96'
        endif
        ny =  96
        if ( .not. present(mask_gg) ) then
          !--- The user has not supplied the name of a file containing
          !--- the land mask. Define one here
          mask_gg_name = "my_landmask_0.25_t95_192_96"
        endif
        !--- Set the output file name here
        if (olaplon) then
          if ( curr_mask_type == -1 ) then
            grid_file_out = 'grid_desc_t95_o_nomask.nc'
          else
            grid_file_out = 'grid_desc_t95_o_'//trim(mask_gg_name)//'.nc'
          endif
        else
          if ( curr_mask_type == -1 ) then
            grid_file_out = 'grid_desc_t95_nomask.nc'
          else
            grid_file_out = 'grid_desc_t95_'//trim(mask_gg_name)//'.nc'
          endif
        endif

      case ( 't127', 'T127' ) !--- T127
        if (olaplon) then
          nx =  257
          grid_name = 'T127 Gaussian Grid -- 257x128'
        else
          nx =  256
          grid_name = 'T127 Gaussian Grid -- 256x128'
        endif
        ny = 128
        if ( .not. present(mask_gg) ) then
          !--- The user has not supplied the name of a file containing
          !--- the land mask. Define one here
          mask_gg_name = "my_landmask_0.25_t127_256_128"
        endif
        !--- Set the output file name here
        if (olaplon) then
          if ( curr_mask_type == -1 ) then
            grid_file_out = 'grid_desc_t127_o_nomask.nc'
          else
            grid_file_out = 'grid_desc_t127_o_'//trim(mask_gg_name)//'.nc'
          endif
        else
          if ( curr_mask_type == -1 ) then
            grid_file_out = 'grid_desc_t127_nomask.nc'
          else
            grid_file_out = 'grid_desc_t127_'//trim(mask_gg_name)//'.nc'
          endif
        endif

      case default
        write(6,'(3a)')'create_grid_file_gg: Invalid grid type --> ',trim(grid_alias),' <--'
        call xit("create_grid_file_gg",-2)
    end select

    if ( len_trim(mask_gg_name).le.0 ) then
      write(6,*)'create_grid_file_gg: A land mask file is required.'
      call xit("create_grid_file_gg",-3)
    endif

    !--- Allocate space for local arrays
    grid_size = nx*ny
    grid_rank = 2
    grid_corners = 4

    allocate( grid_dims(grid_rank) )
    allocate( grid_edges(2,grid_rank) )
    allocate( grid_imask(grid_size) )
    allocate( grid_center_lat(grid_size) )
    allocate( grid_center_lon(grid_size) )
    allocate( grid_corner_lat(grid_corners,grid_size) )
    allocate( grid_corner_lon(grid_corners,grid_size) )
    allocate( grid_area(grid_size) )

    !--- Define the grid dims array written to the grid description file
    grid_dims(1) = nx
    grid_dims(2) = ny

    !-----------------------------------------------------------------------
    !--- Define longitudes of cell centers and cell corners
    !-----------------------------------------------------------------------

    grid_edges(:,:) = 0
    if (olaplon) then
      dlon = 360.0_8/dble(nx-1)
      grid_edges(2,1) = 1
    else
      dlon = 360.0_8/dble(nx)
    endif

    do iy=1,ny
      do ix=1,nx
        atm_add = (iy-1)*nx + ix

        centerlon = dble(ix-1)*dlon
        minlon = centerlon - half*dlon
        maxlon = centerlon + half*dlon

        grid_center_lon(atm_add  ) = centerlon
        !--- These corners must be consistent with the grid_corner_lon array
        grid_corner_lon(1,atm_add) = minlon
        grid_corner_lon(2,atm_add) = maxlon
        grid_corner_lon(3,atm_add) = maxlon
        grid_corner_lon(4,atm_add) = minlon
      end do
    end do

    !--- Determine latitude center and corners using gaussg/trigl

    !--- Determine gaussian latitudes
    call gaussg(ny/2,sinl,wgt,cosl,rad,wocs)
    call trigl (ny/2,sinl,wgt,cosl,rad,wocs)

    !--- Convert rad to degrees
    cfac=180.d0/3.14159265358979323846d0
    glat(1:ny) = rad(1:ny) * cfac

    !--- Determine latitude boundaries
    arithmetic_mean_gg_bnds = .false.
    geometric_mean_gg_bnds  = .true.
    if ( arithmetic_mean_gg_bnds ) then
      !--- Note: tan( (A+B)/2 ) = ( cos(B)-cos(A) )/( sin(A)-sin(B) )
      !--- Therefore the bounds are arithmetic means of adjacent gaussian latitudes
      do iy=2,ny
        glat_bnds(iy) = atan((cos(rad(iy-1)) - cos(rad(iy  )))/ &
                             (sin(rad(iy  )) - sin(rad(iy-1))))
      enddo
    else if ( geometric_mean_gg_bnds ) then
      !--- Put boundaries at the geometric mean of adjacent gaussian latitudes
      do iy=2,ny/2
        glat_bnds(iy) = -sqrt( rad(iy-1)*rad(iy) )
        glat_bnds(ny-iy+2) = -glat_bnds(iy)
      enddo
      glat_bnds(ny/2+1) = 0.0_8
    else
      write(6,'(3a)')'create_grid_file_gg: Invalid gg boundary type.'
      call xit("create_grid_file_gg",-4)
    endif
    !--- Convert to degrees
    glat_bnds(2:ny) = glat_bnds(2:ny) * cfac
    glat_bnds(1)    = -90.0_8
    glat_bnds(ny+1) =  90.0_8

    if ( verbose > 0 ) then
      write(6,'(a)')'Latitude from gaussg/trigl'
      write(6,'(5g24.16)')glat(1:ny)
      write(6,'(a)')'Latitude bounds'
      write(6,'(5g24.16)')glat_bnds(1:ny+1)
      stop
    endif

    do iy=1,ny
      do ix=1,nx
        atm_add = (iy-1)*nx + ix
        grid_center_lat(atm_add  ) = glat(iy)
        !--- These corners must be consistent with the grid_corner_lat array
        grid_corner_lat(1,atm_add) = glat_bnds(iy)
        grid_corner_lat(2,atm_add) = glat_bnds(iy)
        grid_corner_lat(3,atm_add) = glat_bnds(iy+1)
        grid_corner_lat(4,atm_add) = glat_bnds(iy+1)
      enddo
    enddo

!xxx    !-----------------------------------------------------------------------
!xxx    !--- Calculate grid cell areas
!xxx    !-----------------------------------------------------------------------
!xxx    !--- Assign gauss weights normalized st the total area is 1
!xxx    total_area = 0.0_8
!xxx    do iy=1,ny
!xxx      do ix=1,nx
!xxx        idx1 = (iy-1)*nx + ix
!xxx        grid_area(idx1) = wgt(iy)
!xxx        total_area = total_area + wgt(iy)
!xxx      enddo
!xxx    enddo
!xxx
!xxx    do idx1=1,nx*ny
!xxx      grid_area(idx1) = grid_area(idx1)/total_area
!xxx    enddo
!xxx
!xxx    write(6,*)'Total area in gg is ',sum(grid_area(1:nx*ny))

    !-----------------------------------------------------------------------
    !--- define mask
    !-----------------------------------------------------------------------

    !--- grid_imask = 1 means the grid cell will be used in the remapping
    !--- grid_imask = 0 means the grid cell will not be used in the remapping
    grid_imask(:) = 1

    if ( len_trim(mask_gg_name) > 0 ) then
      !--- Do not supply a mask unless one is provided
      call cc_open(1,trim(mask_gg_name),remote,1)

      call getfld2(1,wrk,-1,-1,nc4to8("  GC"),-1,ibuf,maxx,ok)
      if ( .not. ok ) then
        !--- Rewind and look for MASK
        rewind(1)
        call getfld2(1,wrk,-1,-1,nc4to8("MASK"),-1,ibuf,maxx,ok)
      endif

      call cc_close(1,trim(mask_gg_name),remote,1)

      if (.not.ok) then
        write(6,*)'Error reading land mask file'
        call xit("create_grid_file_gg",-5)
      endif

      if ( .not.olaplon ) then
        !--- Remove the overlap longitude from the CCCma mask data
        idx2=0
        do iy=1,ibuf(6)
          do ix=1,ibuf(5)
            idx1 = (iy-1)*ibuf(5) + ix
            if ( mod(idx1,ibuf(5)).eq.0 ) then
              !--- This is the overlap longitude ...ignore it
              cycle
            endif
            idx2=idx2+1
            wrk(idx2) = wrk(idx1)
          end do
        end do

        if (idx2 .ne. grid_size) then
          write(6,*)'Grid size mismatch. grid_size,idx2: ', &
            grid_size,idx2
          call xit("create_grid_file_gg",-6)
        endif
      endif

      if ( curr_mask_type == 2 .or. curr_mask_type == 0 ) then
        !--- mask ocean points with a zero imask value
        where (wrk(1:grid_size).gt.-0.5) grid_imask(1:grid_size) = 0
        write(6,*)'Ocean points are masked.'
      else if ( curr_mask_type == 1 ) then
        !--- mask land points with a zero imask value
        where (wrk(1:grid_size).lt.-0.5) grid_imask(1:grid_size) = 0
        write(6,*)'Land points are masked.'
      endif

    endif

    if ( curr_mask_type == -1 ) then
      !--- Do not mask any grid points regardless of whether or not
      !--- a mask file was supplied
      grid_imask(:) = 1
      mask_gg_name = "NoMask"
    endif

    if ( .false. ) then
      !-----------------------------------------------------------------------
      !--- Store this grid description in memory
      !-----------------------------------------------------------------------
      n_grid_desc_list = n_grid_desc_list + 1
      if ( n_grid_desc_list > max_grid_desc ) then
        write(6,*)'create_grid_file_gg: Exceeded max_grid_desc = ',max_grid_desc
        call xit("create_grid_file_gg",-7)
      endif
      call assign_one_grid( grid_desc_list(n_grid_desc_list) )
      curr_grid = n_grid_desc_list
    endif

    !-----------------------------------------------------------------------
    !--- Create a file containing the current grid description
    !-----------------------------------------------------------------------
    call write_grid_desc_file( grid_file_out, auxfile=trim(mask_gg_name) )

    !-----------------------------------------------------------------------
    !--- Clean up
    !-----------------------------------------------------------------------
    deallocate( grid_dims )
    deallocate( grid_edges )
    deallocate( grid_imask )
    deallocate( grid_center_lat )
    deallocate( grid_center_lon )
    deallocate( grid_corner_lat )
    deallocate( grid_corner_lon )
    deallocate( grid_area )

  end subroutine create_grid_file_gg

  !***********************************************************************
  !--- Create a grid description file for the CCCma ocean grid
  !***********************************************************************
  subroutine create_grid_file_ocn(grid_alias, grid_file_out, curr_grid, olaplon, ocn_desc)

    character(*), intent(in)           :: grid_alias
    character(*), intent(inout)        :: grid_file_out
    integer,      intent(out)          :: curr_grid
    logical,      intent(in)           :: olaplon
    character(*), intent(in), optional :: ocn_desc

    !--- Local
    integer :: ix, iy, idx1, idx2
    integer :: nCell
    integer :: verbose=1
    character(128) :: ocn_datadesc = " "
    character(512) :: strng
    logical :: remote
    logical :: vel_grid

    integer :: ocn_nx, ocn_ny, ocn_nz, ocn_nx_test
    real (ws_r8), allocatable :: ocn_lat(:), ocn_lath(:), ocn_lon(:), ocn_lonh(:)
    real, allocatable :: ocn_z(:), ocn_zh(:), ocn_zbot(:), ocn_beta_t(:), ocn_beta_v(:)

    real (ws_r8) :: cval1, cval2
    real (ws_r8) ::   &
      dLon, dLat,        &! grid spacing in lat, lon
      minLon, maxLon,    &! lon range for each cell
      minLat, maxLat,    &! lat range for each cell
      midLat, midLon      ! latitude,longitude of cell center

    !-----------------------------------------------------------------------
    ! Start executable statements
    !-----------------------------------------------------------------------

    !--- vel_grid = T implies creation of an ocean velocity grid
    !--- vel_grid = F implies creation of an ocean tracer grid
    vel_grid=.false.

    !--- ocn_datadesc is the name of a CCCma format file contianing the ocean grid
    !--- description. Specifically this file must contain the variables
    !---    LON, LONH, LAT, LATH, Z, ZH, ZBOT
    ocn_datadesc = ' '

    !--- Initialize the descriptive grid name that is written to the grid
    !--- description file created below
    grid_name = " "

    select case ( trim(adjustl(grid_alias)) )
      case ( 'ocnt', 'ocnt_256x192' ) !--- Ocean tracer grid at 256x192
        if (olaplon) then
          nx = 257
          grid_name = 'CCCma Ocean Tracer Grid -- 257x192'
        else
          nx = 256
          grid_name = 'CCCma Ocean Tracer Grid -- 256x192'
        endif
        ny =  192
        vel_grid=.false.
        if ( .not. present(ocn_desc) ) then
          !--- Define the name of a file containing the ocean data description
          !--- This file must include: LON, LONH, LAT, LATH, Z, ZH, ZBOT
          ocn_datadesc = "uwl_canesm4_1_1_256x192l40_datadesc"
        endif
        !--- Set the output file name here
        if (olaplon) then
          grid_file_out = 'grid_desc_ocnt_o_'//trim(ocn_datadesc)//'.nc'
        else
          grid_file_out = 'grid_desc_ocnt_'//trim(ocn_datadesc)//'.nc'
        endif

      case ( 'ocnv', 'ocnv_256x192' ) !--- Ocean velocity grid at 256x192
        if (olaplon) then
          nx = 257
          grid_name = 'CCCma Ocean Velocity Grid -- 257x192'
        else
          nx = 256
          grid_name = 'CCCma Ocean Velocity Grid -- 256x192'
        endif
        ny =  192
        vel_grid=.true.
        if ( .not. present(ocn_desc) ) then
          !--- Define the name of a file containing the ocean data description
          !--- This file must include: LON, LONH, LAT, LATH, Z, ZH, ZBOT
          ocn_datadesc = "uwl_canesm4_1_1_256x192l40_datadesc"
        endif
        !--- Set the output file name here
        if (olaplon) then
          grid_file_out = 'grid_desc_ocnv_o_'//trim(ocn_datadesc)//'.nc'
        else
          grid_file_out = 'grid_desc_ocnv_'//trim(ocn_datadesc)//'.nc'
        endif

      case default
        if ( present(ocn_desc) ) then
          !--- Use the supplied ocean data description file and assume a T-grid
          !--- This ocean data desc file will be processed below
          vel_grid=.false.
        else
          write(6,'(3a)')'create_grid_file_ocn: Invalid grid type --> ',trim(grid_alias),' <--'
          call xit("create_grid_file_ocn",-1)
        endif
    end select

    if ( present(ocn_desc) ) then
      !--- Use the ocean data description file provided in the subroutine parameter list
      if ( len_trim(ocn_desc)>0 ) then
        ocn_datadesc = trim(ocn_desc)
        !--- Set the output file name here
        if (vel_grid) then
          if (olaplon) then
            grid_file_out = 'grid_desc_ocnv_o_'//trim(adjustl(ocn_desc))//'.nc'
          else
            grid_file_out = 'grid_desc_ocnv_'//trim(adjustl(ocn_desc))//'.nc'
          endif
        else
          if (olaplon) then
            grid_file_out = 'grid_desc_ocnt_o_'//trim(adjustl(ocn_desc))//'.nc'
          else
            grid_file_out = 'grid_desc_ocnt_'//trim(adjustl(ocn_desc))//'.nc'
          endif
        endif

        !--- Read the grid dimensions from the CCCma ocean data description file
        call cc_open(1, trim(ocn_datadesc), remote, 1)

        !--- Read LAT from ocean data description file
        rewind(1)
        call getfld2(1,wrk,-1,-1,nc4to8(" LAT"),-1,ibuf,maxx,ok)
        if (.not.ok) then
          write(6,*)'Error reading ocean LAT'
          call xit("create_grid_file_ocn",-2)
        endif

        !--- Close the file but leave any attached file in the local dir
        call cc_close(1, trim(ocn_datadesc), .false., 1)

        !--- Set nx and ny
        if (olaplon) then
          nx = ibuf(5)
        else
          nx = ibuf(5)-1
        endif
        ny = ibuf(6)

        if ( len_trim(grid_name) .le. 0 ) then
          write(grid_name,'(2(a,i3))')'CCCma Ocean Tracer Grid -- ',nx,'x',ny
        endif

      else
        write(6,*)'create_grid_file_ocn: The ocean data description file name is blank.'
        call xit("create_grid_file_ocn",-2)
      endif
    endif

    if ( len_trim(ocn_datadesc).le.0 ) then
      write(6,*)'create_grid_file_ocn: An ocean data description file is required.'
      call xit("create_grid_file_ocn",-3)
    endif

    !--- Allocate space for local arrays
    grid_size = nx*ny
    grid_rank = 2
    grid_corners = 4

    allocate( grid_dims(grid_rank) )
    allocate( grid_edges(2,grid_rank) )
    allocate( grid_imask(grid_size) )
    allocate( grid_center_lat(grid_size) )
    allocate( grid_center_lon(grid_size) )
    allocate( grid_corner_lat(grid_corners,grid_size) )
    allocate( grid_corner_lon(grid_corners,grid_size) )
    allocate( grid_area(grid_size) )

    !-----------------------------------------------------------------------
    !  compute longitudes and latitudes of cell centers and corners.
    !-----------------------------------------------------------------------

    grid_dims(1) = nx
    grid_dims(2) = ny

    grid_edges(:,:) = 0
    if (olaplon) then
      ocn_nx_test = nx
      grid_edges(2,1) = 1
    else
      ocn_nx_test = nx+1
    endif

    call cc_open(1,trim(ocn_datadesc),remote,1)

    !--- Read LAT from ocean data description file
    rewind(1)
    call getfld2(1,wrk,-1,-1,nc4to8(" LAT"),-1,ibuf,maxx,ok)
    if (.not.ok) then
      write(6,*)'Error reading ocean LAT'
      call xit("create_grid_file_ocn",-2)
    endif

    !--- After reading the first 2d record set ocn_nx and ocn_ny
    !--- and allocate 2d fields
    ocn_nx=ibuf(5)
    ocn_ny=ibuf(6)
    if (ocn_nx.ne.ocn_nx_test .or. ocn_ny.ne.ny) then
      write(6,*)'Grid size mismatch for ocean LAT.'
      write(6,*)' ibuf(5),ibuf(6): ',ibuf(5),ibuf(6)
      write(6,*)' ocn_nx ,ocn_ny : ',ocn_nx,ocn_ny
      call xit("create_grid_file_ocn",-2)
    endif
    allocate( ocn_lat(ocn_nx*ocn_ny), ocn_lath(ocn_nx*ocn_ny), &
              ocn_lon(ocn_nx*ocn_ny), ocn_lonh(ocn_nx*ocn_ny), &
              ocn_zbot(ocn_nx*ocn_ny) )
    allocate( ocn_beta_t(ocn_nx*ocn_ny), ocn_beta_v(ocn_nx*ocn_ny) )

    !--- Assign the first field
    ocn_lat(1:ocn_nx*ocn_ny) = wrk(1:ocn_nx*ocn_ny)

    !--- Read LATH from ocean data description file
    rewind(1)
    call getfld2(1,wrk,-1,-1,nc4to8("LATH"),-1,ibuf,maxx,ok)
    if (.not.ok) then
      write(6,*)'Error reading ocean LATH'
      call xit("create_grid_file_ocn",-2)
    endif
    if (ibuf(5).ne.ocn_nx .or. ibuf(6).ne.ocn_ny) then
      write(6,*)'Grid size mismatch for ocean LATH.'
      write(6,*)' ibuf(5),ibuf(6): ',ibuf(5),ibuf(6)
      write(6,*)' ocn_nx ,ocn_ny : ',ocn_nx,ocn_ny
      call xit("create_grid_file_ocn",-2)
    endif
    ocn_lath(1:ocn_nx*ocn_ny) = wrk(1:ocn_nx*ocn_ny)

    !--- Read LON from ocean data description file
    rewind(1)
    call getfld2(1,wrk,-1,-1,nc4to8(" LON"),-1,ibuf,maxx,ok)
    if (.not.ok) then
      write(6,*)'Error reading ocean LON'
      call xit("create_grid_file_ocn",-2)
    endif
    if (ibuf(5).ne.ocn_nx .or. ibuf(6).ne.ocn_ny) then
      write(6,*)'Grid size mismatch for ocean LON.'
      write(6,*)' ibuf(5),ibuf(6): ',ibuf(5),ibuf(6)
      write(6,*)' ocn_nx ,ocn_ny : ',ocn_nx,ocn_ny
      call xit("create_grid_file_ocn",-2)
    endif
    ocn_lon(1:ocn_nx*ocn_ny) = wrk(1:ocn_nx*ocn_ny)

    !--- Read LONH from ocean data description file
    rewind(1)
    call getfld2(1,wrk,-1,-1,nc4to8("LONH"),-1,ibuf,maxx,ok)
    if (.not.ok) then
      write(6,*)'Error reading ocean LONH'
      call xit("create_grid_file_ocn",-2)
    endif
    if (ibuf(5).ne.ocn_nx .or. ibuf(6).ne.ocn_ny) then
      write(6,*)'Grid size mismatch for ocean LONH.'
      write(6,*)' ibuf(5),ibuf(6): ',ibuf(5),ibuf(6)
      write(6,*)' ocn_nx ,ocn_ny : ',ocn_nx,ocn_ny
      call xit("create_grid_file_ocn",-2)
    endif
    ocn_lonh(1:ocn_nx*ocn_ny) = wrk(1:ocn_nx*ocn_ny)

    !--- Read Z from ocean data description file
    rewind(1)
    call getfld2(1,wrk,-1,-1,nc4to8("   Z"),-1,ibuf,maxx,ok)
    if (.not.ok) then
      write(6,*)'Error reading ocean Z'
      call xit("create_grid_file_ocn",-2)
    endif

    !--- After reading the first 1d record set ocn_nz and allocate 1d fields
    ocn_nz = ibuf(5)
    allocate( ocn_z(ocn_nz), ocn_zh(ocn_nz) )
    ocn_z(1:ocn_nz) = wrk(1:ocn_nz)

    !--- Read ZH from ocean data description file
    rewind(1)
    call getfld2(1,wrk,-1,-1,nc4to8("  ZH"),-1,ibuf,maxx,ok)
    if (.not.ok) then
      write(6,*)'Error reading ocean ZH'
      call xit("create_grid_file_ocn",-2)
    endif
    if (ibuf(5).ne.ocn_nz) then
      write(6,*)'Grid size mismatch for ocean ZH.'
      write(6,*)' ibuf(5),ibuf(6): ',ibuf(5),ibuf(6)
      write(6,*)' ocn_nz : ',ocn_nz
      call xit("create_grid_file_ocn",-2)
    endif
    ocn_zh(1:ocn_nz) = wrk(1:ocn_nz)

    !--- Read ZBOT from ocean data description file
    rewind(1)
    call getfld2(1,wrk,-1,-1,nc4to8("ZBOT"),-1,ibuf,maxx,ok)
    if (.not.ok) then
      write(6,*)'Error reading ocean ZBOT'
      call xit("create_grid_file_ocn",-2)
    endif
    if (ibuf(5).ne.ocn_nx .or. ibuf(6).ne.ocn_ny) then
      write(6,*)'Grid size mismatch for ocean ZBOT.'
      write(6,*)' ibuf(5),ibuf(6): ',ibuf(5),ibuf(6)
      write(6,*)' ocn_nx ,ocn_ny : ',ocn_nx,ocn_ny
      call xit("create_grid_file_ocn",-2)
    endif
    ocn_zbot(1:ocn_nx*ocn_ny) = wrk(1:ocn_nx*ocn_ny)

    call cc_close(1,trim(ocn_datadesc),remote,1)

    !--- Determine ocean surface beta mask on tracer grid (1 for ocean, 0 for land)
    call betainf(ocn_nx,ocn_ny,0.0,ocn_zh(1),ocn_zbot,ocn_beta_t)

    !--- Determine ocean surface beta mask on velocity grid (1 for ocean, 0 for land)
    call ocn_betau(ocn_beta_t,ocn_beta_v,ocn_nx,ocn_ny)

    if ( .not. olaplon ) then
      !--- Remove the overlap longitude from the 2d CCCma fields
      idx2=0
      do iy=1,ocn_ny
        do ix=1,ocn_nx
          idx1 = (iy-1)*ocn_nx + ix
          if (mod(idx1,ocn_nx).eq.0) then
            !--- This is the overlap longitude ...ignore it
            cycle
          endif
          idx2=idx2+1
          ocn_lat (idx2) = ocn_lat (idx1)
          ocn_lath(idx2) = ocn_lath(idx1)
          ocn_lon (idx2) = ocn_lon (idx1)
          ocn_lonh(idx2) = ocn_lonh(idx1)
          ocn_zbot(idx2) = ocn_zbot(idx1)
          ocn_beta_t(idx2) = ocn_beta_t(idx1)
          ocn_beta_v(idx2) = ocn_beta_v(idx1)
        end do
      end do

      if (idx2 .ne. grid_size) then
        write(6,*)'Grid size mismatch. grid_size,idx2: ',grid_size,idx2
        call xit("create_grid_file_ocn",-2)
      endif
    endif

    !--- Assign lat/lon values using info read from the data description file
    !---
    !--- ocn_lon  represents grid centers for the ocean tracer grid
    !--- ocn_lonh represents grid centers for the ocean velocity grid
    !--- Velocity grid centers are tracer grid cell boundaries and visa versa.
    !--- The ocean tracer grid at the south pole starts with the southern
    !--- edge of the first cell at -90 degrees
    !--- The ocean velocity grid starts with its cell center at northern
    !--- edge of the first tracer grid cell.
    !---
    do iy=1,ny

       !--- lat/lath have the slowest varying dimension so the first nx elements
       !--- of the ocn_lat/ocn_lath all contain the first latitude values.
       !--- These lat values increment every (iy-1)*nx + 1 elements.

       if (vel_grid) then
         !--- Define an ocean velocity grid
         idx1 = (iy-1)*nx + 1
         if (iy.eq.ny) then
           !--- This is the grid cell closest to the north pole
           maxLat = 90.0_ws_r8
           minLat = ocn_lat (idx1)
           ! midLat = ocn_lath(idx1)
           midLat = 0.5_ws_r8*(minLat+maxLat)
         else
           !--- Otherwise use the previous lat value for the southern edge
           idx2 = iy*nx + 1
           maxLat = ocn_lat(idx2)
           midLat = ocn_lath(idx1)
           minLat = ocn_lat (idx1)
         endif
       else
         !--- Define an ocean tracer grid
         idx1 = (iy-1)*nx + 1
         if (iy.eq.1) then
           !--- This is the grid cell closest to the south pole
           minLat = -90.0_ws_r8
         else
           !--- Otherwise use the previous lath value for the southern edge
           idx2 = (iy-2)*nx + 1
           minLat = ocn_lath(idx2)
         endif
         maxLat = ocn_lath(idx1)
         midLat = ocn_lat (idx1)
       endif

       do ix=1,nx
         !--- lon/lonh have the fastest varying dimension so the ocn_lon/ocn_lonh
         !--- values increment in sequence but repeat every nx elements

         if (vel_grid) then
           !--- Define an ocean velocity grid
           idx1 = (iy-1)*nx + ix
           if (ix.eq.nx) then
             !--- This is the grid cell nearest to the last longitude
             idx2 = (iy-1)*nx + 1
             maxLon = ocn_lon(idx2)
           else
             idx2 = (iy-1)*nx + ix + 1
             maxLon = ocn_lon(idx2)
           endif
           minLon = ocn_lon (idx1)
           midLon = ocn_lonh(idx1)
         else
           !--- Define an ocean tracer grid
           idx1 = (iy-1)*nx + ix
           if (ix.eq.1) then
             !--- This is the grid cell nearest to the first longitude
             ! minLon = ocn_lonh(grid_size)
             minLon = ocn_lon(1) - 0.5_ws_r8*(ocn_lon(2)-ocn_lon(1))
             maxLon = ocn_lonh(idx1)
           else if (ix.eq.nx) then
             !--- This is the grid cell nearest to the last longitude
             idx2 = (iy-1)*nx + ix - 1
             minLon = ocn_lonh(idx2)
!xxx             cval1 = ocn_lon(idx1)
!xxx             cval2 = ocn_lon(idx1+1)
!xxx             if (cval1 .le. 0.0_ws_r8) then
!xxx               cval1 = 360.0_ws_r8 + cval1
!xxx             endif
!xxx             if (cval2 .le. 0.0_ws_r8) then
!xxx               cval2 = 360.0_ws_r8 + cval2
!xxx             endif
             maxLon = ocn_lonh(1+idx2)
           else
             idx2 = (iy-1)*nx + ix - 1
             minLon = ocn_lonh(idx2)
             maxLon = ocn_lonh(idx1)
           endif
           midLon = ocn_lon (idx1)
           if (ix.eq.nx) then
             !--- This is the grid cell nearest to the last longitude
             !--- Ensure that the min/max values are not negative
             if (minLon .le. 0.0_ws_r8) then
               minLon = 360.0_ws_r8 + minLon
             endif
             if (maxLon .le. 0.0_ws_r8) then
               maxLon = 360.0_ws_r8 + maxLon
             endif
           endif
         endif

         !--- Convert negative longitudes to longitudes in the range 0-360
         ! if (minLon.lt.0.0_ws_r8 .and. minLon.ge.-360.0_ws_r8) then
         !   minLon = 360.0_ws_r8 + minLon
         ! endif
         ! if (maxLon.lt.0.0_ws_r8 .and. maxLon.ge.-360.0_ws_r8) then
         !   maxLon = 360.0_ws_r8 + maxLon
         ! endif
         ! if (midLon.lt.0.0_ws_r8 .and. midLon.ge.-360.0_ws_r8) then
         !   midLon = 360.0_ws_r8 + midLon
         ! endif

         nCell = (iy-1)*nx + ix

         if ( .false. ) then
           !--- Debugging output ---
           if ( ix == 1 ) then
             !--- This is the grid cell nearest to the first longitude
             write(6,*)'OCN First lon ...ix,iy,nCell: ',ix,iy,nCell
             write(6,*)'OCN First lon ...minLon,midLon,maxLon: ',minLon,midLon,maxLon
             write(6,*)'OCN First lon ...minLat,midLat,maxLat: ',minLat,midLat,maxLat
           endif
           if (ix == nx) then
             !--- This is the grid cell nearest to the last longitude
             write(6,*)'OCN Last lon ...ix,iy,nCell: ',ix,iy,nCell
             write(6,*)'OCN Last lon ...minLon,midLon,maxLon: ',minLon,midLon,maxLon
             write(6,*)'OCN Last lon ...minLat,midLat,maxLat: ',minLat,midLat,maxLat
           endif
         endif

         grid_center_lat(nCell)   = midLat
         grid_corner_lat(1,nCell) = minLat
         grid_corner_lat(2,nCell) = minLat
         grid_corner_lat(3,nCell) = maxLat
         grid_corner_lat(4,nCell) = maxLat

         grid_center_lon(nCell)   = midLon
         grid_corner_lon(1,nCell) = minLon
         grid_corner_lon(2,nCell) = maxLon
         grid_corner_lon(3,nCell) = maxLon
         grid_corner_lon(4,nCell) = minLon
       end do

    end do

    !-----------------------------------------------------------------------
    !--- Calculate grid cell areas
    !-----------------------------------------------------------------------
!xxx    call compute_grid_cell_area(grid_area, grid_corner_lat, grid_corner_lon, gridname="ocn")

    !-----------------------------------------------------------------------
    !--- define mask
    !-----------------------------------------------------------------------

    !--- grid_imask = 1 means the grid cell will be used in the remapping
    !--- grid_imask = 0 means the grid cell will not be used in the remapping
    grid_imask(:) = 1

    if ( curr_mask_type == 2 ) then
      !--- mask ocean points with a zero imask value
      write(6,*)'Ocean points are masked.'
      if (vel_grid) then
        !--- Use ocean velocity grid mask
        where (ocn_beta_v(1:grid_size) > 0.5) grid_imask = 0
      else
        !--- Use ocean tracer grid mask
        where (ocn_beta_t(1:grid_size) > 0.5) grid_imask = 0
      endif
    else if ( curr_mask_type == 1 .or. curr_mask_type == 0 ) then
      !--- mask land points with a zero imask value
      write(6,*)'Land points are masked.'
      if (vel_grid) then
        !--- Use ocean velocity grid mask
        where (ocn_beta_v(1:grid_size) < 0.5) grid_imask = 0
      else
        !--- Use ocean tracer grid mask
        where (ocn_beta_t(1:grid_size) < 0.5) grid_imask = 0
      endif
    else if ( curr_mask_type == -1 ) then
      !--- Do not mask any grid points
      grid_imask = 1
      idx1 = index(grid_file_out,".nc")
      if ( idx1 > 0 ) then
        grid_file_out = grid_file_out(1:idx1-1) // "_nomask.nc"
      endif
      ocn_datadesc = trim(ocn_datadesc) // " NoMask"
    endif

    if (vel_grid) then
      !--- Always mask the northern most velocity grid cell
      do ix=1,nx
        idx1 = (ny-1)*nx + ix
        grid_imask(idx1) = 0
      enddo
    endif

    if ( .false. ) then
      !-----------------------------------------------------------------------
      !--- Store this grid description in memory
      !-----------------------------------------------------------------------
      n_grid_desc_list = n_grid_desc_list + 1
      if ( n_grid_desc_list > max_grid_desc ) then
        write(6,*)'create_grid_file_ocn: Exceeded max_grid_desc = ',max_grid_desc
        call xit("create_grid_file_ocn",-6)
      endif
      call assign_one_grid( grid_desc_list(n_grid_desc_list) )
      curr_grid = n_grid_desc_list
    endif

    !-----------------------------------------------------------------------
    !--- Create a file contianing the current grid description
    !-----------------------------------------------------------------------
    call write_grid_desc_file( grid_file_out, auxfile=trim(ocn_datadesc) )

    if ( verbose > 10 ) then
      !--- Print grid cell centers and corners
      call print_grid()
    endif

    !-----------------------------------------------------------------------
    !  Clean up temporary arrays
    !-----------------------------------------------------------------------
    deallocate( ocn_lat, ocn_lath, ocn_lon, ocn_lonh, &
                 ocn_z, ocn_zh, ocn_zbot, ocn_beta_t, ocn_beta_v )
    deallocate( grid_dims )
    deallocate( grid_edges )
    deallocate( grid_imask )
    deallocate( grid_center_lat )
    deallocate( grid_center_lon )
    deallocate( grid_corner_lat )
    deallocate( grid_corner_lon )
    deallocate( grid_area )

  end subroutine create_grid_file_ocn

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
! usage:  create_nemo_orca input_namelist
!
! input_namelist may define:
!   nemo_coord_file_in  ... input netcdf file
!                           grid description used by NEMO defining
!                             glamt, gphit, tmask :: T grid coords and land mask
!                             glamu, gphiu, umask :: U grid coords and land mask
!                             glamv, gphiv, vmask :: V grid coords and land mask
!                             glamf, gphif, fmask :: F grid coords and land mask
!   nemo_scrip_file_out ... output netcdf files
!                           grid description used by SCRIP defining
!                             grid_imask(grid_size)
!                             grid_center_lat(grid_size)
!                             grid_center_lon(grid_size)
!                             grid_corner_lat(grid_size, grid_corners)
!                             grid_corner_lon(grid_size, grid_corners)
!   nemo_lon            ... name of netcdf variable containing longitudes (e.g. glamt)
!   nemo_lat            ... name of netcdf variable containing latitudes (e.g. gphit)
!   nemo_mask_file      ... name of netcdf file containing the mask variable
!                           if the mask variable is not found in nemo_coord_file_in
!   nemo_mask_vname     ... name of the netcdf mask variable in nemo_mask_file
!   grid_name           ... a string used for the "title" meta data written to
!                           nemo_scrip_file_out
!
! Examples of the input namelist:
!
! 2 degree ORCA2 grid with land masked
! &grid_inputs
!     nemo_coord_file_in  = 'orca2_mesh_mask.nc'
!     nemo_scrip_file_out = 'scrip_input_grid_nemo_orca2.nc'
!     grid_name = 'NEMO ORCA2 T grid -- 182x149'
!     nemo_mask_vname = 'tmask'
!     nemo_lon = 'glamt'
!     nemo_lat = 'gphit'
! /
!
! 2 degree ORCA2 grid without land masked
! &grid_inputs
!     nemo_coord_file_in  = 'orca2_mesh_mask.nc'
!     nemo_scrip_file_out = 'scrip_input_grid_nemo_orca2_nomask.nc'
!     grid_name = 'NEMO ORCA2 T grid -- 182x149 nomask'
!     nemo_mask_vname = 'none'
!     nemo_lon = 'glamt'
!     nemo_lat = 'gphit'
! /
!
! 1 degree ORCA1 grid with land masked
! &grid_inputs
!     nemo_coord_file_in  = 'orca1_mesh_mask.nc'
!     nemo_scrip_file_out = 'scrip_input_grid_nemo_orca1.nc'
!     grid_name = 'NEMO ORCA1 T grid -- 362x292'
!     nemo_mask_vname = 'tmask'
!     nemo_lon = 'glamt'
!     nemo_lat = 'gphit'
! /
!
! 1 degree ORCA1 grid without land masked
! &grid_inputs
!     nemo_coord_file_in  = 'orca1_mesh_mask.nc'
!     nemo_scrip_file_out = 'scrip_input_grid_nemo_orca1.nc'
!     grid_name = 'NEMO ORCA1 T grid -- 362x292 nomask'
!     nemo_mask_vname = 'none'
!     nemo_lon = 'glamt'
!     nemo_lat = 'gphit'
! /
!
! 1/4 degree ORCA025 grid with land masked
! &grid_inputs
!     nemo_coord_file_in  = 'orca025_mesh_mask.nc'
!     nemo_scrip_file_out = 'scrip_input_grid_nemo_orca025.nc'
!     grid_name = 'NEMO ORCA025 T grid -- 1442x1021'
!     nemo_mask_vname = 'tmask'
!     nemo_lon = 'glamt'
!     nemo_lat = 'gphit'
! /
!
! 1/4 degree ORCA025 grid without land masked
! &grid_inputs
!     nemo_coord_file_in  = 'orca025_mesh_mask.nc'
!     nemo_scrip_file_out = 'scrip_input_grid_nemo_orca025_nomask.nc'
!     grid_name = 'NEMO ORCA025 T grid -- 1442x1021 nomask'
!     nemo_mask_vname = 'none'
!     nemo_lon = 'glamt'
!     nemo_lat = 'gphit'
! /


!xxx module create_nemo_orca
!xxx 
!xxx   use netcdf
!xxx 
!xxx   use SCRIP_KindsMod
!xxx   use constants
!xxx   use SCRIP_IOUnitsMod
!xxx   use SCRIP_NetcdfMod
!xxx   use SCRIP_ErrorMod
!xxx 
!xxx   IMPLICIT none
!xxx 
!xxx   !-----------------------------------------------------------------------
!xxx   !     module variables that describe the grid
!xxx 
!xxx   INTEGER (kind=SCRIP_i4), parameter :: &
!xxx        grid_rank = 2,         &
!xxx        grid_corners = 4
!xxx   INTEGER (kind=SCRIP_i4) :: nx, ny, grid_size
!xxx   INTEGER (kind=SCRIP_i4), dimension(2) :: &
!xxx        grid_dims,    &  ! size of x, y dimensions
!xxx        grid_dim_ids     ! ids of the x, y dimensions
!xxx   INTEGER (kind=SCRIP_i4), ALLOCATABLE, DIMENSION(:) :: &
!xxx        grid_imask          ! land-sea mask
!xxx   REAL (kind=SCRIP_i4), ALLOCATABLE, DIMENSION(:) :: &
!xxx        grid_center_lat, &  ! lat/lon coordinates for
!xxx        grid_center_lon     ! each grid center in degrees
!xxx   REAL (kind=SCRIP_r8), ALLOCATABLE, DIMENSION(:,:) :: &
!xxx        grid_corner_lat, &  ! lat/lon coordinates for
!xxx        grid_corner_lon     ! each grid corner in degrees
!xxx   REAL (kind=SCRIP_r8), ALLOCATABLE, DIMENSION(:,:,:) :: &
!xxx        corner_lon, &
!xxx        corner_lat
!xxx   REAL (kind=SCRIP_r8), PARAMETER :: circle = 360.0
!xxx 
!xxx   INTEGER (kind=SCRIP_i4), dimension(4) :: start, count
!xxx 
!xxx   !-----------------------------------------------------------------------
!xxx   !     module variables that describe the netcdf file
!xxx 
!xxx   INTEGER (kind=SCRIP_i4) :: &
!xxx        ncstat,           &   ! general netCDF status variable
!xxx        ncid_in
!xxx 
!xxx   !-----------------------------------------------------------------------
!xxx   !     module variables that describe the grid
!xxx 
!xxx   CHARACTER(SCRIP_charLength) :: grid_name_out = 'Remapped NEMO grid for SCRIP'
!xxx 
!xxx   !-----------------------------------------------------------------------
!xxx   !--- land/ocean mask related variables
!xxx   REAL (kind=SCRIP_r4), ALLOCATABLE, DIMENSION(:,:) :: land_mask
!xxx   INTEGER (kind=SCRIP_i4) :: ncid_mask, varid_mask, ndims_mask, nx_mask, ny_mask
!xxx 
!xxx contains
!xxx 
!xxx   ! ==============================================================================
!xxx   
!xxx   SUBROUTINE convert(nm_in)
!xxx   
!xxx     ! -----------------------------------------------------------------------------
!xxx     ! - input variables
!xxx   
!xxx     CHARACTER(SCRIP_charLength), INTENT(in) ::  &
!xxx       nm_in
!xxx   
!xxx     ! -----------------------------------------------------------------------------
!xxx     ! - local variables
!xxx 
!xxx     CHARACTER(SCRIP_charLength) :: nemo_coord_file_in, nemo_scrip_file_out, &
!xxx         nemo_lon, nemo_lat, corn_lon, corn_lat, grid_name, &
!xxx         nemo_mask_file, nemo_mask_vname
!xxx     INTEGER (kind=SCRIP_i4), dimension(2) :: offset
!xxx     INTEGER (kind=SCRIP_i4) :: iunit
!xxx 
!xxx     namelist /grid_inputs/ nemo_coord_file_in, nemo_scrip_file_out, &
!xxx       nemo_lon, nemo_lat, nemo_mask_file, nemo_mask_vname, grid_name
!xxx   
!xxx     !-----------------------------------------------------------------------
!xxx     ! - namelist describing the processing
!xxx     !   note that mask_value is the minimum good value,
!xxx     !   so that where the mask is less than the value is masked
!xxx 
!xxx     nemo_coord_file_in = "orca1_mesh_mask.nc"
!xxx     nemo_lon = "glamt"
!xxx     nemo_lat = "gphit"
!xxx     nemo_mask_file  = "none"
!xxx     nemo_mask_vname = "tmask"
!xxx     nemo_scrip_file_out = 'scrip_input_grid_nemo_orca1.nc'
!xxx     grid_name = "NEMO ORCA1 T grid -- 362x292"
!xxx 
!xxx     !--- Read the namelist file
!xxx     call SCRIP_IOUnitsGet(iunit)
!xxx     open(iunit, file=nm_in, status='old', form='formatted')
!xxx     read(iunit, nml=grid_inputs) 
!xxx     call SCRIP_IOUnitsRelease(iunit)
!xxx 
!xxx     if (nemo_lon(1:4) .ne. 'glam' .or. nemo_lat(1:4) .ne. 'gphi') then
!xxx       write(6,*) 'lon name does not start with "glam" or lat name does not start with "gphi"'
!xxx       stop
!xxx     endif
!xxx 
!xxx     if ( trim(adjustl(nemo_mask_vname)) == "none" ) then
!xxx       grid_name_out = 'NEMO ORCA1 T grid -- 362x292 nomask'
!xxx     else
!xxx       grid_name_out = 'NEMO ORCA1 T grid -- 362x292'
!xxx     endif
!xxx     if (len_trim(grid_name) > 0) grid_name_out = grid_name
!xxx 
!xxx     ! set up the names of the corner variables for a given input
!xxx     ! the offset represents what needs to be added to (i,j) to get to the correct 
!xxx     ! element in the corner arrays to correspond to the point northeast of the center
!xxx     if (nemo_lon(5:5) == "t") then
!xxx       corn_lon = "glamf"
!xxx       corn_lat = "gphif"
!xxx       offset = (/ 0,0 /)
!xxx     else if (nemo_lon(5:5) == "u") then
!xxx       corn_lon = "glamv"
!xxx       corn_lat = "gphiv"
!xxx       offset = (/ 1,0 /)
!xxx     else if (nemo_lon(5:5) == "v") then
!xxx       corn_lon = "glamu"
!xxx       corn_lat = "gphiu"
!xxx       offset = (/ 0,1 /)
!xxx     else
!xxx       write(6,*) 'unknown nemo_lon name'
!xxx       stop
!xxx     endif
!xxx 
!xxx     write(6,'(a)') "Processing " // trim(nemo_coord_file_in)
!xxx     call convertNEMO(nemo_coord_file_in, nemo_lon, nemo_lat, corn_lon, corn_lat, &
!xxx                      offset, nemo_mask_file, nemo_mask_vname, nemo_scrip_file_out)
!xxx 
!xxx   END SUBROUTINE convert

  !***********************************************************************
  !--- Create a grid description file for a NEMO ORCA grid
  !***********************************************************************
  subroutine create_grid_file_orca(grid_alias, grid_file_out, curr_grid, olaplon, nemo_coords)

    character(*), intent(in)           :: grid_alias
    character(*), intent(inout)        :: grid_file_out
    integer,      intent(out)          :: curr_grid
    logical,      intent(in)           :: olaplon

    !--- nemo_coords is the name of an existing NEMO netcdf file containing ORCA
    !--- coordinate values and land mask in variables named as follows
    !---     glamt, gphit, tmask :: T grid coords and land mask
    !---     glamu, gphiu, umask :: U grid coords and land mask
    !---     glamv, gphiv, vmask :: V grid coords and land mask
    !---     glamf, gphif, fmask :: F grid coords and land mask
    !--- The glam* variables are longitudes
    !--- The gphi* variables are latitudes
    character(*), intent(in), optional :: nemo_coords

    !--- Local
    integer :: ix, iy, idx1, idx2
    integer :: verbose=1
    character(128) :: orca_coords = " "
    character(512) :: strng
    logical :: remote_orca_coords
    character(1) :: nemo_grid_type
    integer(ws_i4) :: offset(2)
    character(32) :: mask_vname
    character(32) :: cent_lon, cent_lat, corn_lon, corn_lat

    !-----------------------------------------------------------------------
    ! Start executable statements
    !-----------------------------------------------------------------------

    nemo_grid_type = 't'

    !--- Define names of the corner and mask variables consistent with cent_(lon|lat)
    !--- The offset represents what needs to be added to (i,j) to get to the correct 
    !--- element in the corner arrays to correspond to the point northeast of the center
    select case ( trim(adjustl(nemo_grid_type)) )
      case ('t', 'T')
        mask_vname = 'tmask'
        cent_lon   = 'glamt'
        cent_lat   = 'gphit'
        corn_lon   = "glamf"
        corn_lat   = "gphif"
        offset     = (/ 0,0 /)
      case ('u', 'U')
        mask_vname = 'umask'
        cent_lon   = 'glamu'
        cent_lat   = 'gphiu'
        corn_lon   = "glamv"
        corn_lat   = "gphiv"
        offset     = (/ 1,0 /)
      case ('v', 'V')
        mask_vname = 'vmask'
        cent_lon   = 'glamv'
        cent_lat   = 'gphiv'
        corn_lon   = "glamu"
        corn_lat   = "gphiu"
        offset     = (/ 0,1 /)
      case default
        write(6,*)'create_grid_file_orca: Invalid NEMO grid type name ',trim(nemo_grid_type)
        call xit("create_grid_file_orca",-1)
    end select

    !--- Initialize the descriptive grid name that is written to the grid
    !--- description file created below
    grid_name = " "
    orca_coords = " "

    select case ( trim(adjustl(grid_alias)) )
      case ( 'orca2', 'ORCA2' ) !--- ORCA2 grid
        !--- grid name is used in the title attribute
        grid_name = 'ORCA2 T grid -- 182x149'

        !--- Define the name of a file containing the ORCA2 coordinates
        orca_coords = "nemo_3.4_orca2_mesh_mask.nc"

        !--- Set the output file name
        grid_file_out = 'grid_desc_orca2_'//trim(orca_coords)

      case ( 'orca1', 'ORCA1' ) !--- ORCA1 grid
        !--- grid name is used in the title attribute
        grid_name = 'ORCA1 T grid -- 362x292'

        !--- Define the name of a file containing the ORCA1 coordinates
        orca_coords = "nemo_3.4_orca1_coordinates_ukorca1_with_mask.nc"

        !--- Set the output file name here
        grid_file_out = 'grid_desc_orca1_'//trim(orca_coords)

      case ( 'orca025', 'ORCA025' ) !--- ORCA025 grid
        !--- grid name is used in the title attribute
        grid_name = 'ORCA025 T grid -- 1442x1021'

        !--- Define the name of a file containing the ORCA1 coordinates
        orca_coords = "nemo3.4.1_orca025_mesh_mask_ngb.nc"

        !--- Set the output file name here
        grid_file_out = 'grid_desc_orca025_'//trim(orca_coords)

      case default
        if ( .not. present(nemo_coords) ) then
          write(6,'(3a)')'create_grid_file_orca: Invalid grid type --> ',trim(grid_alias),' <--'
          call xit("create_grid_file_orca",-1)
        endif

    end select

    if ( present(nemo_coords) ) then
      if ( len_trim(nemo_coords)>0 ) then
        !--- Use the NEMO coordinates file provided in the subroutine parameter list
        orca_coords = trim(nemo_coords)

        !--- grid name is used in the title attribute
        grid_name = 'ORCA grid'

        !--- Set the output file name
        grid_file_out = 'grid_desc_orca_'//trim(adjustl(nemo_coords))
      else
        write(6,*)'create_grid_file_orca: The grid description file name is blank.'
        call xit("create_grid_file_orca",-2)
      endif
    endif

    if ( len_trim(orca_coords).le.0 ) then
      write(6,*)'create_grid_file_orca: A file containing ORCA coordinate values is required.'
      call xit("create_grid_file_orca",-3)
    endif

    !--- The orca_coords file must exist in the cwd before convertnemo is called
    if ( .not. locate_file(trim(orca_coords),remote_orca_coords) ) then
      write(6,*)'create_grid_file_orca: Missing coordinates file ',trim(orca_coords)
      call xit("create_grid_file_orca",-4)
    endif

    !-----------------------------------------------------------------------
    !  Compute longitudes and latitudes of cell centers and corners.
    !-----------------------------------------------------------------------

    call convertNEMO(orca_coords, cent_lon, cent_lat, corn_lon, corn_lat, &
                         offset, mask_vname, grid_file_out)

!xxx    !--- Compute grid cell area
!xxx    call compute_grid_cell_area( grid_area, grid_corner_lat, grid_corner_lon, &
!xxx                                gridname=trim(grid_name) )

    if ( verbose > 10 ) then
      !--- Print grid cell centers and corners
      call print_grid()
    endif

    if ( curr_mask_type == -1 ) then
      !--- No grid points are masked
      idx1 = index(grid_file_out,".nc")
      if ( idx1 > 0 ) then
        grid_file_out = grid_file_out(1:idx1-1) // "_nomask.nc"
      endif
      orca_coords = trim(orca_coords)//"  NoMask"
    endif

    !-----------------------------------------------------------------------
    !--- Create a file containing the current grid description
    !-----------------------------------------------------------------------
    call write_grid_desc_file( grid_file_out, auxfile=trim(orca_coords) )

    if (remote_orca_coords) then
      strng=' '
      strng="release "//trim(orca_coords)
      call system(trim(strng))
    endif

    !-----------------------------------------------------------------------
    !  Clean up temporary arrays
    !-----------------------------------------------------------------------
    deallocate( grid_dims )
    deallocate( grid_edges )
    deallocate( grid_imask )
    deallocate( grid_center_lat )
    deallocate( grid_center_lon )
    deallocate( grid_corner_lat )
    deallocate( grid_corner_lon )
    deallocate( grid_area )
    if ( allocated(nav_lon) ) deallocate(nav_lon)
    if ( allocated(nav_lat) ) deallocate(nav_lat)

  end subroutine create_grid_file_orca

  ! ==============================================================================

  subroutine convertNEMO(grid_file_in, cent_lon, cent_lat, corn_lon, corn_lat, &
                         offset, nemo_mask_vname, grid_file_out)

    !--- Convert a NEMO coordinates.nc file to a SCRIP format grid file.

    !--- Names of coordinate variables found in the input netcdf file
    character(*), intent(in) :: cent_lon, cent_lat, corn_lon, corn_lat

    !--- offset is what needs to be added to (i,j) to get to the correct element
    !--- in the corner arrays to correspond to the point northeast of the center
    integer (kind=ws_i4), intent(in), dimension(2) :: offset

    !--- The name of the input netcdf file containing NEMO format coords
    character(*), intent(in) :: grid_file_in

    !--- The name of the output netcdf file containing SCRIP format coords
    character(*), intent(in) :: grid_file_out

    !--- The name of the mask variable found in the input netcdf file
    character(*), intent(in) :: nemo_mask_vname

    !--- grid coordinates and masks

    !--- center longitude
    real(ws_r8), allocatable, dimension(:,:) :: glam
    !--- center latitude
    real(ws_r8), allocatable, dimension(:,:) :: gphi
    !--- corner longitude
    real(ws_r8), allocatable, dimension(:,:) :: glamc
    !--- corner latitude
    real(ws_r8), allocatable, dimension(:,:) :: gphic
    !--- expanded corner arrays
    real(ws_r8), allocatable, dimension(:,:) :: expand_clon, expand_clat

    real(ws_r8), allocatable :: gwrk1(:)

    integer(ws_i4) :: grid_dim_ids(2)
    real(ws_r8), allocatable, dimension(:,:,:) :: corner_lon, corner_lat
    real(ws_r8), parameter :: circle = 360.0_ws_r8
    integer(ws_i4) :: nc_start(4), nc_count(4)
    integer(ws_i4) :: ncstat, ncid_in
    real(ws_r4), allocatable :: land_mask(:,:)
    integer(ws_i4), allocatable :: iland_mask(:,:)
    integer(ws_i4) :: ncid_mask, varid_mask, ndims_mask, nx_mask, ny_mask
    integer(ws_i4) :: varid_nav_lon, varid_nav_lat

    !--- other local variables

    integer (kind=ws_i4) :: i, j, n, iunit, im1, jm1, imid, isame, ic, jc
    integer (kind=ws_i4) :: varid_lam, varid_phi, varid_lamc, varid_phic
    integer (kind=ws_i4) :: jdim
    logical :: found
    integer (kind=ws_i4), dimension(4) :: grid_dimids  ! input fields have 4 dims
    real (kind=ws_r8) :: tmplon, dxt, dyt
    character(len=128) :: afmt, bfmt

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

    !--- For NEMO input grids, assume that variable names are glam, glamc etc.
    !--- Assume that 1st 2 dimensions of these variables are x and y directions.
    !--- These assumptions are made by NEMO, so should be valid for coordinates.nc.

    ncstat = nf90_open( grid_file_in, NF90_NOWRITE, ncid_in )
    call handle_nc_err(ncstat, ncid_in, msg="open nemo grid file")

    !--- find dimids for 'glam'
    !--- use dimids to get dimlengths
    !--- allocate glam array
    !--- get glam from file

    !--- Allocate the array to hold nx,ny dimension sizes
    grid_rank = 2
    grid_corners = 4
    allocate( grid_dims(grid_rank) )
    allocate( grid_edges(2,grid_rank) )
    grid_edges(:,:) = 0
    grid_edges(2,1) = 2

    !--- Determine netcdf variable ids for input coordinate variable names
    ncstat = nf90_inq_varid( ncid_in, trim(cent_lon), varid_lam )
    call handle_nc_err(ncstat, ncid_in)
    ncstat = nf90_inq_varid( ncid_in, trim(corn_lon), varid_lamc )
    call handle_nc_err(ncstat, ncid_in)
    ncstat = nf90_inq_varid( ncid_in, trim(cent_lat), varid_phi )
    call handle_nc_err(ncstat, ncid_in)
    ncstat = nf90_inq_varid( ncid_in, trim(corn_lat), varid_phic )
    call handle_nc_err(ncstat, ncid_in)

    ncstat = nf90_inq_varid( ncid_in, "nav_lon", varid_nav_lon )
    call handle_nc_err(ncstat, ncid_in)
    ncstat = nf90_inq_varid( ncid_in, "nav_lat", varid_nav_lat )
    call handle_nc_err(ncstat, ncid_in)

    !--- Assign grid_dims assuming that the first 2 dimensions of any coordinate
    !--- variable will be X and Y (X varying the fastest)
    ncstat = nf90_inquire_variable( ncid_in, varid_lam, dimids=grid_dimids(:) )
    call handle_nc_err(ncstat, ncid_in)
    do jdim = 1,2
       ncstat = nf90_inquire_dimension( ncid_in, grid_dimids(jdim), len=grid_dims(jdim) )
       call handle_nc_err(ncstat, ncid_in)
    end do
    write(6,'("input grid dimensions are:",2i6)') grid_dims

    !--- Allocate space for local arrays
    nx = grid_dims(1)
    ny = grid_dims(2)
    grid_size = nx*ny
    allocate( grid_imask(grid_size) )
    allocate( grid_center_lat(grid_size) )
    allocate( grid_center_lon(grid_size) )
    allocate( grid_corner_lat(grid_corners,grid_size) )
    allocate( grid_corner_lon(grid_corners,grid_size) )
    allocate( grid_area(grid_size) )
    grid_area(:) = 0.0_ws_r8

    !--- Assume that all coordinate dimensions are the same and allocate
    allocate( glam(nx,ny), glamc(nx,ny), gphi(nx,ny), gphic(nx,ny) )

    !--- Read coordinate values from the input netcdf file
    ncstat = nf90_get_var( ncid_in, varid_lam,  glam(:,:) )
    call handle_nc_err(ncstat, ncid_in)
    ncstat = nf90_get_var( ncid_in, varid_lamc, glamc(:,:) )
    call handle_nc_err(ncstat, ncid_in)
    ncstat = nf90_get_var( ncid_in, varid_phi,  gphi(:,:) )
    call handle_nc_err(ncstat, ncid_in)
    ncstat = nf90_get_var( ncid_in, varid_phic, gphic(:,:) )
    call handle_nc_err(ncstat, ncid_in)

    allocate( nav_lon(nx,ny), nav_lat(nx,ny) )
    ncstat = nf90_get_var( ncid_in, varid_nav_lon, nav_lon(:,:) )
    call handle_nc_err(ncstat, ncid_in)
    ncstat = nf90_get_var( ncid_in, varid_nav_lat, nav_lat(:,:) )
    call handle_nc_err(ncstat, ncid_in)

    if (.false.) then
      !--- Write representative center and vertex coordinates to a text file
      open(97,file='nemo_edge_coords.txt',form='formatted')
      afmt = " "
      bfmt = '3(3x,"(",f10.4,",",f10.4,")")'
      write(afmt,'(5a)')'(a,i5,'//trim(bfmt)//'," ... ",'//trim(bfmt)//')'
      do j=1,ny
        write(97,trim(afmt)) "cent: i=1,2,3, ... nx-2,nx-1,nx j=",j, &
          glam(1,j),gphi(1,j),glam(2,j),gphi(2,j),glam(3,j),gphi(3,j), &
          glam(nx-2,j),gphi(nx-2,j), glam(nx-1,j),gphi(nx-1,j), glam(nx,j),gphi(nx,j)
        write(97,trim(afmt)) "vert: i=1,2,3, ... nx-2,nx-1,nx j=",j, &
          glamc(1,j),gphic(1,j),glamc(2,j),gphic(2,j),glamc(3,j),gphic(3,j), &
          glamc(nx-2,j),gphic(nx-2,j), glamc(nx-1,j),gphic(nx-1,j), glamc(nx,j),gphic(nx,j)
        write(97,*)" "
      enddo

      do i=1,nx
        write(97,trim(afmt)) "cent: j=1,2,3, ... ny-2,ny-1,ny i=",i, &
          glam(i,1),gphi(i,1),glam(i,2),gphi(i,2),glam(i,3),gphi(i,3), &
          glam(i,ny-2),gphi(i,ny-2), glam(i,ny-1),gphi(i,ny-1), glam(i,ny),gphi(i,ny)
        write(97,trim(afmt)) "vert: j=1,2,3, ... ny-2,ny-1,ny i=",i, &
          glamc(i,1),gphic(i,1),glamc(i,2),gphic(i,2),glamc(i,3),gphic(i,3), &
          glamc(i,ny-2),gphic(i,ny-2), glamc(i,ny-1),gphic(i,ny-1), glamc(i,ny),gphic(i,ny)
        write(97,*)" "
      enddo
      close(97)
      write(6,*)'Created nemo_edge_coords.txt'
    endif
    
    !-----------------------------------------------------------------------
    !--- Define the land/ocean mask (land=0, ocean=1)
  
    !--- initialize mask as all ocean
    grid_imask(:) = 1

    !--- Look for the mask in the input nemo coord description file
    ncstat = nf90_inq_varid( ncid_in, trim(nemo_mask_vname), varid_mask )
    if (ncstat /= NF90_NOERR) then
      !--- If the mask variable does not exist then issue a warning and move on
      write(6,*)'convertNEMO: mask variable ',trim(nemo_mask_vname), &
                ' not found in input file ',trim(grid_file_in)
      write(6,*)'convertNEMO: Output mask will include all grid points.'

    else
      !--- Otherwise read the mask from the input file
      ncstat = nf90_inquire_variable(ncid_in, varid_mask, ndims=ndims_mask)
      call handle_nc_err(ncstat, ncid_in, msg="inquire nemo mask variable Ndims")
      grid_dimids(:) = -1
      ncstat = nf90_inquire_variable(ncid_in, varid_mask, dimids=grid_dimids(1:ndims_mask) )
      call handle_nc_err(ncstat, ncid_in, msg="inquire nemo mask variable dimIDs")

      !--- Assume that the first 2 dimensions are nx,ny
      ncstat = nf90_inquire_dimension( ncid_in, grid_dimids(1), len=nx_mask )
      call handle_nc_err(ncstat, ncid_in)
      ncstat = nf90_inquire_dimension( ncid_in, grid_dimids(2), len=ny_mask )
      call handle_nc_err(ncstat, ncid_in)
      if ( nx /= nx_mask .or. ny /= ny_mask ) then
        write(6,*)'convertNEMO: Land_mask dimensions do not match NEMO grid'
        write(6,'(a,i4,a,i4)')'  nx_mask = ',nx_mask,'  ny_mask = ',ny_mask
        write(6,'(a,i4,a,i4)')'  nx_nemo = ',nx,'  ny_nemo = ',ny
        stop
      endif

      !--- Define the portion of the netcdf mask variable to be read below
      nc_start(:) = 1
      nc_count(:) = (/ nx_mask, ny_mask, one_i4, one_i4 /)

      !--- Read land_mask from the user supplied file
      allocate( land_mask(nx_mask,ny_mask) )
      ncstat = nf90_get_var( ncid_in, varid_mask, land_mask(:,:), start=nc_start, count=nc_count )
      call handle_nc_err(ncstat, ncid_in)

      if ( .false. ) then
        !--- Print the grid mask (0 = cell not used, 1 = cell used)
        allocate( iland_mask(nx_mask,ny_mask) )
        iland_mask = nint(land_mask)
        do j=ny_mask,1,-1
          write(6,'(i4,2x,500i1)')j,iland_mask(:,j)
        enddo
        call flush(6)
        deallocate( iland_mask )
      endif

      !--- Assign the 1D grid_imask array using the mask read from the file
      n=0
      do j=1,ny_mask
        do i=1,nx_mask
          n=n+1
          if ( curr_mask_type == 2 ) then
            !--- mask ocean points assuming land_mask == 0 over land
            if (land_mask(i,j) /= 0) then
              grid_imask(n) = 0
            else
              grid_imask(n) = 1
            endif
          else
            !--- Use the mask that was read from the file
            grid_imask(n) = land_mask(i,j)
          endif
        enddo
      enddo

      deallocate( land_mask )

    endif

    if ( curr_mask_type == -1 ) then
      !--- Do not mask any grid points even if a mask variable
      !--- exists in the grid description file
      grid_imask(:) = 1
    endif

    !--- Close the nemo coordinates file
    ncstat = nf90_close(ncid_in)
    call handle_nc_err(ncstat, ncid_in)
  
    !-----------------------------------------------------------------------
    ! Corners must be ordered as follows:
    !     4 3
    !     1 2
    !
    ! Create a pair of arrays for the corners where expand_clon(1,1) corresponds
    ! to the south west corner of a box centered at (glam(1,1),gphi(1,1))
  
    allocate ( expand_clon(nx+1,ny+1), expand_clat(nx+1,ny+1) )

    expand_clon(:,:) = spval
    expand_clat(:,:) = spval

    !--- First the internal points
    do j = 2,ny
      do i = 2,nx
        ic = i + offset(1) - 1
        jc = j + offset(2) - 1
        expand_clon(i,j) = glamc(ic,jc)
        expand_clat(i,j) = gphic(ic,jc)
      enddo
    enddo

    !--- Then the boundary points
    !--- Assume that in the data read from the coordinates file there are 2 overlap
    !--- longitudes so that the following conditions hold (glam == lon, gphi == lat)
    !---     glam(1,:) == glam(nx-1,:)     gphi(1,:) == gphi(nx-1,:)
    !---     glam(2,:) == glam(nx,:)       gphi(2,:) == gphi(nx,:)
    !--- and
    !---     glamc(1,:) == glamc(nx-1,:)   gphic(1,:) == gphic(nx-1,:)
    !---     glamc(2,:) == glamc(nx,:)     gphic(2,:) == gphic(nx,:)
    do i=2,nx
      ic = i + offset(1) - 1
      !--- First (south/bottom) row
      expand_clon(i,1) = glamc(ic,1)
      expand_clat(i,1) = gphic(ic,1) - (gphic(ic,2)-gphic(ic,1))

      !--- Last (north/top) row
      expand_clon(i,ny+1) = glamc(ic,ny)
      expand_clat(i,ny+1) = gphic(ic,ny)
    enddo

    do j=2,ny
      jc = j + offset(2) - 1
      !--- First (west/left) column
      expand_clon(1,j) = glamc(nx-2,jc)
      expand_clat(1,j) = gphic(nx-2,jc)

      !--- Last (east/right) column
      expand_clon(nx+1,j) = glamc(2,jc)
      expand_clat(nx+1,j) = gphic(2,jc)
    enddo

    !--- Then the corner points
    ic = nx-2
    expand_clon(1,1) = glamc(ic,1)
    expand_clat(1,1) = gphic(ic,1) - (gphic(ic,2)-gphic(ic,1))

    ic = 2
    expand_clon(nx+1,1) = glamc(ic,1)
    expand_clat(nx+1,1) = gphic(ic,1) - (gphic(ic,2)-gphic(ic,1))

    imid = (nx-1)/2 + 1
    ic = 2*imid - nx + 1
    expand_clon(1,ny+1) = glamc(ic,ny)
    expand_clat(1,ny+1) = gphic(ic,ny)

    ic = 2*imid
    expand_clon(nx+1,ny+1) = glamc(ic,ny)
    expand_clat(nx+1,ny+1) = gphic(ic,ny)

    if ( .true.) then
      write(6,'(a,2f10.4)')'Bottom left  (1,1)    corner on grid: lon,lat = ', &
                            expand_clon(1,1),expand_clat(1,1)
      write(6,'(a,2f10.4)')'Bottom right (nx,1)   corner on grid: lon,lat = ', &
                            expand_clon(nx+1,1),expand_clat(nx+1,1)
      write(6,'(a,2f10.4)')'Top    left  (1,ny)   corner on grid: lon,lat = ', &
                            expand_clon(1,ny+1),expand_clat(1,ny+1)
      write(6,'(a,2f10.4)')'Top    right (nx,ny)  corner on grid: lon,lat = ', &
                            expand_clon(nx+1,ny+1),expand_clat(nx+1,ny+1)
    endif

    !--- Verify that no corner points were missed
    found=.false.
    if ( any(expand_clon == spval) ) then
      found=.true.
      jdim = count( expand_clon .eq. spval )
      write(6,*)'expand_clon contains ',jdim,' spvals'
      do i=1,nx+1
        do j=1,ny+1
          if (expand_clon(i,j) == spval) write(6,*)i,j,' expand_clon = spval'
        enddo
      enddo
    endif
    if ( any(expand_clat == spval) ) then
      found=.true.
      jdim = count( expand_clat .eq. spval )
      write(6,*)'expand_clat contains ',jdim,' spvals'
      do i=1,nx+1
        do j=1,ny+1
          if (expand_clat(i,j) == spval) write(6,*)i,j,' expand_clat = spval'
        enddo
      enddo
    endif
    if (found) then
      write(6,*)'Missing values were found in the corner array.'
      call xit("convertnemo",-1)
    endif

    !--- Assign corner arrays on grids shaped arrays
    allocate ( corner_lon(4,nx,ny), corner_lat(4,nx,ny) )

    corner_lon(:,:,:) = zero
    corner_lat(:,:,:) = zero

    ! top-right corner
    corner_lon(3,:,:) = expand_clon(2:nx+1,2:ny+1)
    corner_lat(3,:,:) = expand_clat(2:nx+1,2:ny+1)
  
    ! top-left corner
    corner_lon(4,:,:) = expand_clon(1:nx,2:ny+1)
    corner_lat(4,:,:) = expand_clat(1:nx,2:ny+1)
  
    ! bottom-right corner 
    corner_lon(2,:,:) = expand_clon(2:nx+1,1:ny)
    corner_lat(2,:,:) = expand_clat(2:nx+1,1:ny)
  
    ! bottom-left corner
    corner_lon(1,:,:) = expand_clon(1:nx,1:ny)
    corner_lat(1,:,:) = expand_clat(1:nx,1:ny)

    !--- Zero out small absolute values
    do i=1,size(corner_lon,2)
      do j=1,size(corner_lon,3)
        if (abs(corner_lon(1,i,j)) < 1e-20) corner_lon(1,i,j) = zero
        if (abs(corner_lon(2,i,j)) < 1e-20) corner_lon(2,i,j) = zero
        if (abs(corner_lon(3,i,j)) < 1e-20) corner_lon(3,i,j) = zero
        if (abs(corner_lon(4,i,j)) < 1e-20) corner_lon(4,i,j) = zero
      enddo
    enddo
    do i=1,size(corner_lat,2)
      do j=1,size(corner_lat,3)
        if (abs(corner_lat(1,i,j)) < 1e-20) corner_lat(1,i,j) = zero
        if (abs(corner_lat(2,i,j)) < 1e-20) corner_lat(2,i,j) = zero
        if (abs(corner_lat(3,i,j)) < 1e-20) corner_lat(3,i,j) = zero
        if (abs(corner_lat(4,i,j)) < 1e-20) corner_lat(4,i,j) = zero
      enddo
    enddo

!    do i=1,size(corner_lat,2)
!      do j=1,size(corner_lat,3)
!        if (abs(corner_lon(1,i,j)-12.030076) < 1e-5 .and. &
!            abs(corner_lat(1,i,j)-30.011358) < 1e-5) then
!          !--- found (12.030076,  30.011358) ...possible corner of concave element
!          write(*,*)'Concave corner 1 (SW), i,j = ',i,j
!        endif
!        if (abs(corner_lon(2,i,j)-12.030076) < 1e-5 .and. &
!            abs(corner_lat(2,i,j)-30.011358) < 1e-5) then
!          !--- found (12.030076,  30.011358) ...possible corner of concave element
!          write(*,*)'Concave corner 2 (SW), i,j = ',i,j
!        endif
!        if (abs(corner_lon(3,i,j)-12.030076) < 1e-5 .and. &
!            abs(corner_lat(3,i,j)-30.011358) < 1e-5) then
!          !--- found (12.030076,  30.011358) ...possible corner of concave element
!          write(*,*)'Concave corner 3 (SW), i,j = ',i,j
!        endif
!        if (abs(corner_lon(4,i,j)-12.030076) < 1e-5 .and. &
!            abs(corner_lat(4,i,j)-30.011358) < 1e-5) then
!          !--- found (12.030076,  30.011358) ...possible corner of concave element
!          write(*,*)'Concave corner 4 (SW), i,j = ',i,j
!        endif
!      enddo
!    enddo
!
!    do i=156,159
!      do j=94,97
!        call print_one_cell(glam, gphi, corner_lon, corner_lat, i, j)
!      enddo
!    enddo

!    if (.false.) then
!      !--- correct for 0,2pi longitude crossings
!      !--- In practice this means putting all corners into 0,2pi range
!      !--- and ensuring that no box corners are miles from each other.
!      corner_lon(:,:,:) = modulo( corner_lon(:,:,:), circle )
!      do n = 2, grid_corners
!         where    ( corner_lon(n,:,:) - corner_lon(n-1,:,:) < -three*circle*0.25 )
!            corner_lon(n,:,:) = corner_lon(n,:,:) + circle
!         elsewhere( corner_lon(n,:,:) - corner_lon(n-1,:,:) >  three*circle*0.25 )
!            corner_lon(n,:,:) = corner_lon(n,:,:) - circle
!         end where
!      end do
!    endif

    !--- Reshape center and corner arrays for SCRIP format
    grid_center_lon(:) = reshape( glam(:,:), (/ grid_size /) )
    grid_center_lat(:) = reshape( gphi(:,:), (/ grid_size /) )

    deallocate( glam, gphi, glamc, gphic )

    grid_corner_lon(:,:) = reshape( corner_lon(:,:,:), (/ grid_corners, grid_size /) )
    grid_corner_lat(:,:) = reshape( corner_lat(:,:,:), (/ grid_corners, grid_size /) )

    deallocate( corner_lon, corner_lat )

    if ( .false. ) then
      !--- Reassign all corner longitudes to be in the range 0,2pi
      ! grid_center_lon(:)   = modulo( grid_center_lon(:),   circle )
      grid_corner_lon(:,:) = modulo( grid_corner_lon(:,:), circle )

      !--- Ensure that box corners are not too far apart
      do n = 2, grid_corners
         where    ( grid_corner_lon(n,:) - grid_corner_lon(n-1,:) < -three*circle*0.25 )
            grid_corner_lon(n,:) = grid_corner_lon(n,:) + circle
         elsewhere( grid_corner_lon(n,:) - grid_corner_lon(n-1,:) >  three*circle*0.25 )
            grid_corner_lon(n,:) = grid_corner_lon(n,:) - circle
         end where
      end do
    endif

    if ( .false. ) then
      do n=1,nx
        call print_cell_coords(n)
      enddo
      do n=grid_size-nx+1,grid_size
        call print_cell_coords(n)
      enddo
    endif

    contains

      subroutine print_cell_coords(n)
        integer(kind=ws_i4) :: n, i
        write(6,'(i8,2x,"lon,lat = ","(",f11.6,",",f11.6,")  corners =",4(1x,"(",f11.6,",",f11.6,")"))') &
        n,grid_center_lon(n),grid_center_lat(n),(grid_corner_lon(i,n),grid_corner_lat(i,n),i=1,4)
      end subroutine print_cell_coords
  
  end subroutine convertNEMO


  !*****************************************************************************
  !--- Write a previously defined grid description to a netcdf file
  !--- The following variables must be defined prior to this call
  !---    grid_size          ...number of cells in the grid (e.g. nx*ny)
  !---    grid_rank          ...rank of the grid (e.g. 2)
  !---    grid_corners       ...number of grid cell corners (e.g. 4)
  !---    grid_name          ...descriptive name for the grid
  !---    grid_dims(:)       ...shape of the grid (e.g. (nx,ny) )
  !---    grid_imask(:)      ...integer land mask (0 or 1)
  !---    grid_center_lat(:) ...lat coords for each grid center (degrees)
  !---    grid_center_lon(:) ...lon coords for each grid center (degrees)
  !---    grid_corner_lat(:,:) ...lat coords for each grid cell corner (degrees)
  !---    grid_corner_lon(:,:) ...lon coords for each grid cell corner (degrees)
  !---    grid_area(:)       ...area of each grid cell
  !*****************************************************************************
  subroutine write_grid_desc_file(grid_file_out, auxfile)

    character(*), intent(in) :: grid_file_out

    character(*), intent(in), optional :: auxfile

    !--- ncstat            ...general netCDF status variable
    !--- nc_grid_id        ...netCDF grid dataset id
    !--- nc_gridsize_id    ...netCDF grid size dim id
    !--- nc_gridcorn_id    ...netCDF grid corner dim id
    !--- nc_gridrank_id    ...netCDF grid rank dim id
    !--- nc_griddims_id    ...netCDF grid dimension size id
    !--- nc_gridedges_id    ...netCDF grid dimension overlap id
    !--- nc_grdcntrlat_id  ...netCDF grid center lat id
    !--- nc_grdcntrlon_id  ...netCDF grid center lon id
    !--- nc_grdimask_id    ...netCDF grid mask id
    !--- nc_grdcrnrlat_id  ...netCDF grid corner lat id
    !--- nc_grdcrnrlon_id  ...netCDF grid corner lon id
    !--- nc_grdarea_id     ...netCDF grid area id
    integer (kind=4) :: ncstat, &
                        nc_grid_id, &
                        nc_gridsize_id,  &
                        nc_gridcorn_id,  &
                        nc_gridrank_id,  &
                        nc_griddims_id,  &
                        nc_gridedges_id, &
                        nc_edges_id,     &
                        nc_grdcntrlat_id, &
                        nc_grdcntrlon_id, &
                        nc_grdimask_id,   &
                        nc_grdcrnrlat_id, &
                        nc_grdcrnrlon_id, &
                        nc_nav_lon_id,    &
                        nc_nav_lat_id,    &
                        nc_grdarea_id

    !--- netCDF dim id array for 2-d arrays
    integer (kind=4), dimension(2) :: nc_dims2_id
    integer (kind=4), parameter :: edges = 2

    !-----------------------------------------------------------------------
    !--- Open the netCDF file, overwrite any existing file of the same name
    !-----------------------------------------------------------------------

    ncstat = nf90_create (grid_file_out, NF90_CLOBBER, nc_grid_id)
    call handle_nc_err(ncstat, nc_grid_id, msg="call nf90_create...")

    !-----------------------------------------------------------------------
    !--- Set up attributes for netCDF file
    !-----------------------------------------------------------------------

    !***
    !*** Add a descriptive name for the grid in the global meta data "title"
    !***

    ncstat = nf90_put_att (nc_grid_id, NF90_GLOBAL, 'title', trim(grid_name))
    call handle_nc_err(ncstat, nc_grid_id)

    if ( present(auxfile) ) then
      !***
      !*** Include the name of any auxillary data file used
      !***
      ncstat = nf90_put_att (nc_grid_id, NF90_GLOBAL, 'external_input_file', trim(auxfile))
      call handle_nc_err(ncstat, nc_grid_id)
    endif

    !***
    !*** define grid size dimension
    !***

    ncstat = nf90_def_dim (nc_grid_id, 'grid_size', grid_size, nc_gridsize_id)
    call handle_nc_err(ncstat, nc_grid_id)

    !***
    !*** define grid corner dimension
    !***

    ncstat = nf90_def_dim (nc_grid_id, 'grid_corners', grid_corners, nc_gridcorn_id)
    call handle_nc_err(ncstat, nc_grid_id)

    !***
    !*** define grid rank dimension
    !***

    ncstat = nf90_def_dim (nc_grid_id, 'grid_rank', grid_rank, nc_gridrank_id)
    call handle_nc_err(ncstat, nc_grid_id)

    !***
    !*** define grid edge number (always 2)
    !***

    ncstat = nf90_def_dim (nc_grid_id, 'edges', edges, nc_edges_id)
    call handle_nc_err(ncstat, nc_grid_id)

    !***
    !*** define grid dimension size array
    !***

    ncstat = nf90_def_var (nc_grid_id, 'grid_dims', NF90_INT, &
                       nc_gridrank_id, nc_griddims_id)
    call handle_nc_err(ncstat, nc_grid_id)

    !***
    !*** define grid overlap dimension array
    !***

    nc_dims2_id(1) = nc_edges_id
    nc_dims2_id(2) = nc_gridrank_id

    ncstat = nf90_def_var (nc_grid_id, 'grid_edges', NF90_INT, &
                       nc_dims2_id, nc_gridedges_id)
    call handle_nc_err(ncstat, nc_grid_id)

    !***
    !*** define grid center latitude array
    !***

    ncstat = nf90_def_var (nc_grid_id, 'grid_center_lat', NF90_DOUBLE, &
                       nc_gridsize_id, nc_grdcntrlat_id)
    call handle_nc_err(ncstat, nc_grid_id)

    ncstat = nf90_put_att (nc_grid_id, nc_grdcntrlat_id, 'units', 'degrees')
    call handle_nc_err(ncstat, nc_grid_id)

    !***
    !*** define grid center longitude array
    !***

    ncstat = nf90_def_var (nc_grid_id, 'grid_center_lon', NF90_DOUBLE, &
                       nc_gridsize_id, nc_grdcntrlon_id)
    call handle_nc_err(ncstat, nc_grid_id)

    ncstat = nf90_put_att (nc_grid_id, nc_grdcntrlon_id, 'units', 'degrees')
    call handle_nc_err(ncstat, nc_grid_id)

    !***
    !*** define grid mask
    !***

    ncstat = nf90_def_var (nc_grid_id, 'grid_imask', NF90_INT, &
                           nc_gridsize_id, nc_grdimask_id)
    call handle_nc_err(ncstat, nc_grid_id)

    ncstat = nf90_put_att (nc_grid_id, nc_grdimask_id, 'units', 'unitless')
    call handle_nc_err(ncstat, nc_grid_id)

    !***
    !*** define grid corner latitude array
    !***

    nc_dims2_id(1) = nc_gridcorn_id
    nc_dims2_id(2) = nc_gridsize_id

    ncstat = nf90_def_var (nc_grid_id, 'grid_corner_lat', NF90_DOUBLE, &
                           nc_dims2_id, nc_grdcrnrlat_id)
    call handle_nc_err(ncstat, nc_grid_id)

    ncstat = nf90_put_att (nc_grid_id, nc_grdcrnrlat_id, 'units', 'degrees')
    call handle_nc_err(ncstat, nc_grid_id)

    !***
    !*** define grid corner longitude array
    !***

    ncstat = nf90_def_var (nc_grid_id, 'grid_corner_lon', NF90_DOUBLE, &
                           nc_dims2_id, nc_grdcrnrlon_id)
    call handle_nc_err(ncstat, nc_grid_id)

    ncstat = nf90_put_att (nc_grid_id, nc_grdcrnrlon_id, 'units', 'degrees')
    call handle_nc_err(ncstat, nc_grid_id)

    if ( allocated(nav_lon) ) then
      !***
      !*** define nav_lon array
      !***

      ncstat = nf90_def_var (nc_grid_id, 'nav_lon', NF90_FLOAT, &
                         nc_gridsize_id, nc_nav_lon_id)
      call handle_nc_err(ncstat, nc_grid_id)

      ncstat = nf90_put_att (nc_grid_id, nc_nav_lon_id, 'units', 'degrees_east')
      call handle_nc_err(ncstat, nc_grid_id)
    endif

    if ( allocated(nav_lat) ) then
      !***
      !*** define nav_lat array
      !***

      ncstat = nf90_def_var (nc_grid_id, 'nav_lat', NF90_FLOAT, &
                         nc_gridsize_id, nc_nav_lat_id)
      call handle_nc_err(ncstat, nc_grid_id)

      ncstat = nf90_put_att (nc_grid_id, nc_nav_lat_id, 'units', 'degrees_north')
      call handle_nc_err(ncstat, nc_grid_id)
    endif

    if ( allocated(grid_area) .and. write_grid_area ) then
      !***
      !*** define grid area array
      !***

      ncstat = nf90_def_var (nc_grid_id, 'grid_area', NF90_DOUBLE, &
                         nc_gridsize_id, nc_grdarea_id)
      call handle_nc_err(ncstat, nc_grid_id)

      ncstat = nf90_put_att (nc_grid_id, nc_grdarea_id, 'units', 'unitless')
      call handle_nc_err(ncstat, nc_grid_id)

      ncstat = nf90_put_att (nc_grid_id, nc_grdarea_id, 'comment', &
               'Normalized to the unit sphere')
      call handle_nc_err(ncstat, nc_grid_id)
    endif

    !***
    !*** end definition stage
    !***

    ncstat = nf90_enddef(nc_grid_id)
    call handle_nc_err(ncstat, nc_grid_id)

    !-----------------------------------------------------------------------
    !--- Write grid data
    !-----------------------------------------------------------------------

    ncstat = nf90_put_var(nc_grid_id, nc_griddims_id, grid_dims)
    call handle_nc_err(ncstat, nc_grid_id)

    ncstat = nf90_put_var(nc_grid_id, nc_gridedges_id, grid_edges)
    call handle_nc_err(ncstat, nc_grid_id)

    ncstat = nf90_put_var(nc_grid_id, nc_grdimask_id, grid_imask)
    call handle_nc_err(ncstat, nc_grid_id)

    ncstat = nf90_put_var(nc_grid_id, nc_grdcntrlat_id, grid_center_lat)
    call handle_nc_err(ncstat, nc_grid_id)

    ncstat = nf90_put_var(nc_grid_id, nc_grdcntrlon_id, grid_center_lon)
    call handle_nc_err(ncstat, nc_grid_id)

    ncstat = nf90_put_var(nc_grid_id, nc_grdcrnrlat_id, grid_corner_lat)
    call handle_nc_err(ncstat, nc_grid_id)

    ncstat = nf90_put_var(nc_grid_id, nc_grdcrnrlon_id, grid_corner_lon)
    call handle_nc_err(ncstat, nc_grid_id)

    if ( allocated(nav_lon) ) then
      ncstat = nf90_put_var(nc_grid_id, nc_nav_lon_id, &
                            reshape(nav_lon(:,:),(/ grid_size /)) )
      call handle_nc_err(ncstat, nc_grid_id)
    endif

    if ( allocated(nav_lat) ) then
      ncstat = nf90_put_var(nc_grid_id, nc_nav_lat_id, &
                            reshape(nav_lat(:,:),(/ grid_size /)) )
      call handle_nc_err(ncstat, nc_grid_id)
    endif

    if ( allocated(grid_area) .and. write_grid_area ) then
      ncstat = nf90_put_var(nc_grid_id, nc_grdarea_id, grid_area)
      call handle_nc_err(ncstat, nc_grid_id)
    endif

    !-----------------------------------------------------------------------
    !--- Close the netcdf file
    !-----------------------------------------------------------------------

    ncstat = nf90_close(nc_grid_id)
    call handle_nc_err(ncstat, nc_grid_id)

    write(6,'(2a)')'Created grid description file ',trim(grid_file_out)

  end subroutine write_grid_desc_file

  !***********************************************************************
  !--- Read an existing grid description file and store the data
  !--- internally, namely in the array grid_desc_list
  !***********************************************************************
  function read_xxx_grid_desc_file(file_name) result(grid_idx)

    !--- The return value is the index in grid_desc_list pointing to the
    !--- stored data read from the input file
    integer :: grid_idx

    !--- Name of file containing remap data
    character(*), intent(in) :: file_name

    !--- local
    integer (kind=4) :: ncstat, fid, vid
    integer (kind=4), parameter :: edges = 2
    integer :: idx, ic
    logical :: found
    character (len=NF90_MAX_NAME) :: vname
    character (len=128) :: external_input_file

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

    !--- Open the grid description file
    ncstat = nf90_open(file_name, NF90_NOWRITE, fid)
    call handle_nc_err(ncstat, fid)

    !--- Read global meta data
    ncstat = nf90_get_att(fid,NF90_GLOBAL,'title',grid_name)
    call handle_nc_err(ncstat, fid)

    ncstat = nf90_get_att(fid,NF90_GLOBAL,'external_input_file',external_input_file)
    call handle_nc_err(ncstat, fid)

    !--- Read dimension sizes
    ncstat = nf90_inq_dimid(fid, 'grid_size', vid)
    call handle_nc_err(ncstat, fid)
    ncstat = nf90_inquire_dimension(fid, vid, len=grid_size)
    call handle_nc_err(ncstat, fid)

    ncstat = nf90_inq_dimid(fid, 'grid_corners', vid)
    call handle_nc_err(ncstat, fid)
    ncstat = nf90_inquire_dimension(fid, vid, len=grid_corners)
    call handle_nc_err(ncstat, fid)

    ncstat = nf90_inq_dimid(fid, 'grid_rank', vid)
    call handle_nc_err(ncstat, fid)
    ncstat = nf90_inquire_dimension(fid, vid, len=grid_rank)
    call handle_nc_err(ncstat, fid)

    !--- Read grid_dims first so that it may be used in the following
    !--- test to find out if this grid has been read before
    if ( allocated(grid_dims) ) deallocate(grid_dims)
    allocate( grid_dims(grid_rank) )

    !--- Read the shape of the grid into grid_dims
    ncstat = nf90_inq_varid(fid, 'grid_dims', vid)
    call handle_nc_err(ncstat, fid)
    ncstat = nf90_get_var(fid, vid, grid_dims)
    call handle_nc_err(ncstat, fid)

    if ( n_grid_desc_list > 0 ) then
      ic = 0
      do idx=1,n_grid_desc_list
        found = trim(grid_desc_list(idx)%title) .eq. trim(grid_name) .and. &
                grid_desc_list(idx)%size        .eq. grid_size       .and. &
                grid_desc_list(idx)%corners     .eq. grid_corners    .and. &
                grid_desc_list(idx)%rank        .eq. grid_rank
        if ( found ) then
          !--- Also require that all the grid dimensions match
          found = all(grid_desc_list(idx)%dims .eq. grid_dims)
        endif
        if ( found ) then
          !--- This grid description matches that found in the file
          ic = idx
          exit
        endif
      enddo
      if ( ic > 0 ) then
        !--- This file has been seen before
        if ( allocated(grid_dims) ) deallocate(grid_dims)
        grid_idx = ic
        ncstat = nf90_close(fid)
        call handle_nc_err(ncstat, fid)
        return
      endif
    endif

    !--- Read the number of overlap grid points in each dim into grid_edges
    if ( allocated(grid_edges) ) deallocate(grid_edges)
    allocate( grid_edges(edges,grid_rank) )
    ncstat = nf90_inq_varid(fid, 'grid_edges', vid)
    call handle_nc_err(ncstat, fid)
    ncstat = nf90_get_var(fid, vid, grid_edges)
    call handle_nc_err(ncstat, fid)

    !--- Allocate space for variables to be read
    if ( allocated(grid_center_lat) ) deallocate(grid_center_lat)
    if ( allocated(grid_center_lon) ) deallocate(grid_center_lon)
    if ( allocated(grid_corner_lat) ) deallocate(grid_corner_lat)
    if ( allocated(grid_corner_lon) ) deallocate(grid_corner_lon)
    if ( allocated(grid_imask) )      deallocate(grid_imask)
    if ( allocated(grid_area) )       deallocate(grid_area)
    allocate( grid_center_lat(grid_size) )
    allocate( grid_center_lon(grid_size) )
    allocate( grid_corner_lat(grid_corners,grid_size) )
    allocate( grid_corner_lon(grid_corners,grid_size) )
    allocate( grid_imask(grid_size) )
    allocate( grid_area(grid_size) )

    !--- Read the center lat/lon arrays
    ncstat = nf90_inq_varid(fid, 'grid_center_lat', vid)
    call handle_nc_err(ncstat, fid)
    ncstat = nf90_get_var(fid, vid, grid_center_lat)
    call handle_nc_err(ncstat, fid)

    ncstat = nf90_inq_varid(fid, 'grid_center_lon', vid)
    call handle_nc_err(ncstat, fid)
    ncstat = nf90_get_var(fid, vid, grid_center_lon)
    call handle_nc_err(ncstat, fid)

    !--- Read the corner lat/lon arrays
    ncstat = nf90_inq_varid(fid, 'grid_corner_lat', vid)
    call handle_nc_err(ncstat, fid)
    ncstat = nf90_get_var(fid, vid, grid_corner_lat)
    call handle_nc_err(ncstat, fid)

    ncstat = nf90_inq_varid(fid, 'grid_corner_lon', vid)
    call handle_nc_err(ncstat, fid)
    ncstat = nf90_get_var(fid, vid, grid_corner_lon)
    call handle_nc_err(ncstat, fid)

    !--- Read the imask array
    ncstat = nf90_inq_varid(fid, 'grid_imask', vid)
    call handle_nc_err(ncstat, fid)
    ncstat = nf90_get_var(fid, vid, grid_imask)
    call handle_nc_err(ncstat, fid)

    if ( read_grid_area ) then
      !--- Read the area array
      ncstat = nf90_inq_varid(fid, 'grid_area', vid)
      call handle_nc_err(ncstat, fid)
      ncstat = nf90_get_var(fid, vid, grid_area)
      call handle_nc_err(ncstat, fid)
    endif

    !--- Close the netcdf file
    ncstat = nf90_close(fid)
    call handle_nc_err(ncstat, fid)

    !--- Store the name of the file read in grid_disk_file so that it will be
    !--- saved in grid_desc_list along with the rest of the grid info
    grid_disk_file = trim(file_name)

    !--- Store this grid description in memory and assign the array index
    !--- pointing to its location to the function value
    n_grid_desc_list = n_grid_desc_list + 1
    if ( n_grid_desc_list > max_grid_desc ) then
      write(6,*)'read_xxx_grid_desc_file: Exceeded max_grid_desc = ',max_grid_desc
      call xit("read_xxx_grid_desc_file",-1)
    endif
    call assign_one_grid( grid_desc_list(n_grid_desc_list) )
    grid_idx = n_grid_desc_list

    !--- Deallocate all arrays that were used above
    if ( allocated(grid_dims) )       deallocate(grid_dims)
    if ( allocated(grid_center_lat) ) deallocate(grid_center_lat)
    if ( allocated(grid_center_lon) ) deallocate(grid_center_lon)
    if ( allocated(grid_corner_lat) ) deallocate(grid_corner_lat)
    if ( allocated(grid_corner_lon) ) deallocate(grid_corner_lon)
    if ( allocated(grid_imask) )      deallocate(grid_imask)
    if ( allocated(grid_area) )       deallocate(grid_area)
    grid_disk_file = " "

    write(6,'(3a,i8,a)')'read_xxx_grid_desc_file: Read ',trim(file_name), &
        ' into index ',grid_idx,' of grid_desc_list.'

  end function read_xxx_grid_desc_file

  !***********************************************************************
  !--- Handle errors returned from calls to netcdf routines
  !***********************************************************************
  subroutine handle_nc_err(status, fid, msg)

    implicit none

    !--- Input
    integer(kind=4), intent(in) :: status         ! return code from nf90_* call
    integer(kind=4), intent(in) :: fid            ! netCDF file id
    character(*),    intent(in), optional :: msg ! Info string to print on error

    !--- Local
    integer(kind=4) :: ret

    !----- If no error, do nothing
    if (status.eq.NF90_NOERR) return

    if ( present(msg) ) then
      if ( len_trim(msg)>0 ) write(6,*)trim(msg)
    endif
    print *, '*** ',status,' netCDF error: ', nf90_strerror(status)
    ret = nf90_close (fid)

    call xit("handle_nc_err",-1)

  end subroutine handle_nc_err

  !***********************************************************************
  !--- This subroutine finds the roots (in theta) and gaussian weights 
  !--- associated with the legendre polynomial of degree l > 1.
  !***********************************************************************
  subroutine gquad(l,root,w)

    implicit none

    !-----------------------------------------------------------------------
    !
    !     intent(in)
    !
    !-----------------------------------------------------------------------

    integer (kind=ws_i4), intent(in) :: l

    !-----------------------------------------------------------------------
    !
    !     intent(out)
    !
    !-----------------------------------------------------------------------

    real (kind=ws_r8), dimension(l), intent(out) :: root, w

    !-----------------------------------------------------------------------
    !
    !     local
    !
    !-----------------------------------------------------------------------

    integer (kind=ws_i4) :: l1, l2, l22, l3, k, i, j

    real (kind=ws_r8) :: del,co,p1,p2,p3,t1,t2,slope,s,c,pp1,pp2,p00

    !-----------------------------------------------------------------------
    !
    !     Define useful constants.
    !
    !-----------------------------------------------------------------------

    del= pi/float(4*l)
    l1 = l+1
    co = float(2*l+3)/float(l1**2)
    p2 = 1.0
    t2 = -del
    l2 = l/2
    k = 1
    p00 = one/sqrt(two)

    !-----------------------------------------------------------------------
    !
    !     Start search for each root by looking for crossing point.
    !
    !-----------------------------------------------------------------------

    do i=1,l2
   10  t1 = t2
       t2 = t1+del
       p1 = p2
       s = sin(t2)
       c = cos(t2)
       pp1 = 1.0
       p3 = p00
       do j=1,l1
          pp2 = pp1
          pp1 = p3
          p3 = 2.0*sqrt((float(j**2)-0.250)/float(j**2))*c*pp1- &
               sqrt(float((2*j+1)*(j-1)*(j-1))/ &
               float((2*j-3)*j*j))*pp2
       end do
       p2 = pp1
       if ((k*p2).gt.0) goto 10

    !-----------------------------------------------------------------------
    !
    !        Now converge using Newton-Raphson.
    !
    !-----------------------------------------------------------------------

       k = -k
   20  continue
          slope = (t2-t1)/(p2-p1)
          t1 = t2
          t2 = t2-slope*p2
          p1 = p2
          s = sin(t2)
          c = cos(t2)
          pp1 = 1.0
          p3 = p00
          do j=1,l1
             pp2 = pp1
             pp1 = p3
             p3 = 2.0*sqrt((float(j**2)-0.250)/float(j**2))*c*pp1- &
                  sqrt(float((2*j+1)*(j-1)*(j-1))/ &
                  float((2*j-3)*j*j))*pp2
          end do
          p2 = pp1
       if (abs(p2).gt.1.e-10) goto 20
       root(i) = t2
       w(i) = co*(sin(t2)/p3)**2
    end do

    !-----------------------------------------------------------------------
    !
    !     If l is odd, take care of odd point.
    !
    !-----------------------------------------------------------------------

    l22 = 2*l2
    if (l22 .ne. l) then
       l2 = l2+1
       t2 = pi/2.0
       root(l2) = t2
       s = sin(t2)
       c = cos(t2)
       pp1 = 1.0
       p3 = p00
       do j=1,l1
          pp2 = pp1
          pp1 = p3
          p3 = 2.0*sqrt((float(j**2)-0.250)/float(j**2))*c*pp1- &
               sqrt(float((2*j+1)*(j-1)*(j-1))/ &
               float((2*j-3)*j*j))*pp2
       end do
       p2 = pp1
       w(l2) = co/p3**2
    endif

    !-----------------------------------------------------------------------
    !
    !     Use symmetry to compute remaining roots and weights.
    !
    !-----------------------------------------------------------------------

    l3 = l2+1
    do i=l3,l
       root(i) = pi-root(l-i+1)
       w(i) = w(l-i+1)
    end do

  end subroutine gquad

  subroutine compute_grid_cell_area(area, cornerLat_in, cornerLon_in, gridname)
    implicit none

    !--- This routine computes a grid cell area based on corner lat/lon
    !--- coordinates.  It is provided in the case that a user supplied
    !--- area is not available.

    !--- computed area for each grid cell
    real (ws_r8), dimension(:), intent(out) :: area

    !--- latitude  of each cell corner
    real (ws_r8), dimension(:,:), intent(in) :: cornerLat_in

    !--- longitude of each cell corner
    real (ws_r8), dimension(:,:), intent(in) :: cornerLon_in

    !--- Name of the grid (optional)
    character(*), optional :: gridname

    !--- Local
    integer (ws_i4) :: numCells    !--- number of grid cells
    integer (ws_i4) :: numCorners  !--- number of corners in each cell
    integer (ws_i4) :: nCell       !--- loop index for grid cells
    integer (ws_i4) :: nCorner     !--- loop index for corners in each cell
    integer (ws_i4) :: nextCorner  !--- next corner around cell perimeter
    real (ws_r8)    :: dphi        !--- delta (longitude) for this segment

    !--- conversion for deg to rads
    real (ws_r8), parameter :: deg2rad = pi/180.0_ws_r8
    real (ws_r8), dimension(:,:), allocatable :: cornerLon, cornerLat

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

    !---  determine size of grid and initialize
    numCells   = size(CornerLat_in, dim=2)
    numCorners = size(CornerLat_in, dim=1)

    !--- Input corner(lat|lon) arrays are in degrees ...convert to radians
    allocate( cornerLon(numCorners,numCells) )
    allocate( cornerLat(numCorners,numCells) )
    cornerLon(:,:) = deg2rad * CornerLon_in(:,:)
    cornerLat(:,:) = deg2rad * CornerLat_in(:,:)

    !--- Compute area for each cell by integrating around cell edges
    Area(:) = zero
    do nCell=1,numCells

       Area(nCell) = 0.0_ws_r8

       do nCorner=1,numCorners
          nextCorner = mod(nCorner,numCorners) + 1

          !*** trapezoid rule - delta(Lon) is -0.5*dx
          dphi = CornerLon(   nCorner,nCell) -  &
                 CornerLon(nextCorner,nCell)
          if (dphi > pi) then
             dphi = dphi - pi2
          else if (dphi < -pi) then
             dphi = dphi + pi2
          endif
          dphi = 0.5_ws_r8*dphi

          Area(nCell) = Area(nCell) +  &
                        dphi*(sin(CornerLat(   nCorner,nCell)) + &
                              sin(CornerLat(nextCorner,nCell)))
       end do
    end do

    deallocate( cornerLon, cornerLat)

    !--- If the grid cells cover the entire sphere then the sum of all areas
    !--- will be 4*pi (ie the surface area of the unit sphere)
    ! write(6,*)'Normalized total area of grid = ',sum(Area)/(4.0_ws_r8 * pi)
    if ( present(gridname) ) then
      write(6,*)'Total area of grid ',trim(gridname),' = ',sum(Area)
    else
      write(6,*)'Total area of this grid = ',sum(Area)
    endif

  end subroutine compute_grid_cell_area


  !***********************************************************************
  !--- This subroutine takes 1D version of OGCM beta mask (1=ocean,0=land)
  !--- and returns the beta mask for the UV grid.
  !***********************************************************************
  subroutine ocn_betau(beta,betau,imc,jmc)

    implicit none

    real :: beta(*), betau(*)
    integer :: imc,jmc,imcm1,jmcm1,ij,i,j,j1,jm1 

    !--- Modified to use in this module Larry Solheim, May, 2014
    !--- Created WGL July 14, 2010

    !Figure 1: Assumed below is that FORTRAN variable BETA(IJ) =  B(IJ)
    !          and output BETAU(IJ) = BUV(IJ) in this figure, where
    !          IJ = I+(J-1)*IMC
    !
    !          |              |              |              | 
    !      ----+--------------+--------------+--------------+-----------
    !          |              |              |              | 
    !          |              |              |              | 
    !          |              |              |              | 
    !          |              |              |              | 
    !          |              |              |              | 
    !          |              |              |              | 
    !          |                                            | 
    !      ----+--------  BUV(IJ-1) ----  BUV(IJ)  ---------+-----------
    !          |                                            | 
    !          |              |              |              | 
    !          |              |              |              | 
    !          |              |    B(IJ)     |   B(IJ+1)    | 
    !          |              |              |              | 
    !          |              |              |              | 
    !          |                                            | 
    !      ----+-------  BUV(IJ-IMC-1)   BUV(IJ-IMC)   -----+-----------
    !          |                                            | 
    !          |              |              |              | 
    !          |              |              |              | 
    !          |              |  B(IJ-IMC)   | B(IJ-IMC+1)  | 
    !          |              |              |              | 
    !          |              |              |              | 
    !          |              |              |              | 
    !      ----+--------------+--------------+--------------+-----------

    IMCM1 = IMC - 1
    JMCM1 = JMC - 1

    !----   First do all latitudes except for the north pole

    DO ij = 1,IMC*JMC
       betau(IJ) = 2
    ENDDO

    DO J   = 1, JMCM1
       JM1 = J - 1
       DO I = 1, IMCM1
          IJ = JM1*IMC + I
          BETAU(IJ)=BETA(IJ)*BETA(IJ+1)*BETA(IJ+IMC)*BETA(IJ+IMC+1)
       ENDDO
       BETAU(JM1*IMC+IMC) = BETAU(JM1*IMC+1)
    ENDDO

    !----   Ensure all North Pole points are considered ocean

    DO I = 1,IMC
       BETAU(JMCM1*IMC+I) = 1.0
    ENDDO

  end subroutine ocn_betau

end module grid_desc
