module ccc_netcdf
  !***************************************************************
  ! Common netcdf related routines
  !***************************************************************

  use com_cpl
  use grid_comm
  use ncdf

  implicit none

  public

  interface resize_a2d
    module procedure resize_a2d_i4
    module procedure resize_a2d_i8
    module procedure resize_a2d_r4
    module procedure resize_a2d_r8
  end interface

  private :: resize_a2d_i4
  private :: resize_a2d_i8
  private :: resize_a2d_r4
  private :: resize_a2d_r8

  contains

    !***************************************************************************
    !--- Determine the type of grid description file give the file name
    !--- The decision is based on the names of variables in the file
    !***************************************************************************
    function grid_desc_file_type(file_name) result(file_type)
      !--- The name of a grid description file
      character(*), intent(in)  :: file_name

      !--- An integer flag to indicate file type
      integer :: file_type

      !--- Local
      character(64), pointer :: var_names(:)
      integer :: idx1, idx2
      
      !--- Set SCRIP format to be the default grid type
      file_type=1

      !--- Read variable names from the grid description file and use these names
      !--- to identify the type of file that is to be read
      call ncdf_inquire_vars(file_name, var_names=var_names)

      write(6,'(3a,200(2x,a))')'grid_desc_file_type: ',trim(file_name), &
          '   var_names = ',(trim(var_names(idx1)),idx1=1,size(var_names))
      call flush(6)

      idx2=0
      do idx1=1,size(var_names)
        if ( trim(adjustl(var_names(idx1))) .eq. "glamt" ) idx2=idx2+1
        if ( trim(adjustl(var_names(idx1))) .eq. "glamu" ) idx2=idx2+1
        if ( trim(adjustl(var_names(idx1))) .eq. "glamv" ) idx2=idx2+1
        if ( trim(adjustl(var_names(idx1))) .eq. "glamf" ) idx2=idx2+1
        if ( trim(adjustl(var_names(idx1))) .eq. "gphit" ) idx2=idx2+1
        if ( trim(adjustl(var_names(idx1))) .eq. "gphiu" ) idx2=idx2+1
        if ( trim(adjustl(var_names(idx1))) .eq. "gphiv" ) idx2=idx2+1
        if ( trim(adjustl(var_names(idx1))) .eq. "gphif" ) idx2=idx2+1
        if ( trim(adjustl(var_names(idx1))) .eq. "tmask" ) idx2=idx2+1
        if ( trim(adjustl(var_names(idx1))) .eq. "umask" ) idx2=idx2+1
        if ( trim(adjustl(var_names(idx1))) .eq. "vmask" ) idx2=idx2+1
        if ( trim(adjustl(var_names(idx1))) .eq. "fmask" ) idx2=idx2+1
      enddo
      if ( associated(var_names) ) deallocate(var_names)
      if ( idx2 == 8 .or. idx2 == 12 ) then
        !--- This file contains all of the following variables:
        !---     glamt, glamu, glamv, glamf,
        !---     gphit, gphiu, gphiv, gphif
        !---     tmask, umask, vmask, fmask
        !--- Assume it contains NEMO coordinates (e.g. a NEMO mesh_mask file)
        file_type = 2
      endif

      write(6,'(3a,i3)')'grid_desc_file_type: ',trim(file_name), &
                        '   file_type = ',file_type
      call flush(6)

    end function grid_desc_file_type

    !***************************************************************************
    !--- Read specific data from a SCRIP format grid description file
    !***************************************************************************
    subroutine read_scrip_grid_file(file_name, olap, overlap, &
        title, grid_size, grid_dims, grid_edges, imask, use_mask, &
        mid_lon, mid_lat, corn_lon, corn_lat, &
        mid_lon_1d, mid_lat_1d, corn_lon_1d, corn_lat_1d)

      implicit none

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

      !--- Number of overlap longitudes (0,1,2) in the first dimension
      !--- If olap > 0 then the last olap columns will be removed from
      !--- all 2D fields that are read from the grid description file
      !--- and these truncated fields will be returned to the calling
      !--- program as imask, mid_lon, mid_lat, corn_lon, corn_lat
      integer, intent(in), optional :: olap

      !--- If overlap = F then do not include overlaping grid cells in any
      !--- dimension in the output grid description
      logical, intent(in), optional :: overlap

      !--- The grid title read from the grid description file
      character(*), intent(inout), optional :: title

      !--- The total number of points in the grid
      integer, intent(out), optional :: grid_size

      !--- The shape of the grid ie (nx,ny)
      integer, intent(out), optional :: grid_dims(2)

      !--- The number of overlap grids (lower and upper) in each dimension
      integer, intent(out), optional :: grid_edges(2,2)

      !--- Integer mask (1 means vaild grid cell, 0 means ignore this grid cell)
      integer,      intent(out), pointer, optional ::    imask(:,:)

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

      !--- lon/lat of grid cell centers
      real(kind=8), intent(out), pointer, optional ::  mid_lon(:,:)
      real(kind=8), intent(out), pointer, optional ::  mid_lat(:,:)
      !--- Coords as 1D only makes sense for certain grids
      !--- (e.g. regular cartesian, gaussian, etc)
      real(kind=8), intent(out), pointer, optional ::  mid_lon_1d(:)
      real(kind=8), intent(out), pointer, optional ::  mid_lat_1d(:)

      !--- lon/lat of grid cell corners
      !--- It is assumed that there are 4 corners per cell
      real(kind=8), intent(out), pointer, optional :: corn_lon(:,:)
      real(kind=8), intent(out), pointer, optional :: corn_lat(:,:)
      !--- Coords as 1D only makes sense for certain grids
      !--- (e.g. regular cartesian, gaussian, etc)
      real(kind=8), intent(out), pointer, optional :: corn_lon_1d(:)
      real(kind=8), intent(out), pointer, optional :: corn_lat_1d(:)

      !--- Local
      integer :: local_olap
      logical :: keep_overlap
      character(256) :: local_title
      integer(kind=4) :: ncstat, ncfid, ncvid, grid_rank, grid_corners
      integer(kind=4) :: edges
      integer(kind=4) :: local_grid_size
      integer(kind=4) :: local_grid_dims(2), nx, ny
      integer(kind=4), allocatable :: local_grid_edges(:,:)
      integer(kind=4), allocatable :: grid_imask(:)
      real(kind=8), allocatable :: grid_center_lon(:)
      real(kind=8), allocatable :: grid_center_lat(:)
      real(kind=8), allocatable :: grid_corner_lon(:,:)
      real(kind=8), allocatable :: grid_corner_lat(:,:)
      real(kind=8), allocatable :: corn_lon1d(:), corn_lat1d(:)
      real(kind=8), allocatable :: corn_lon_tmp(:,:), corn_lat_tmp(:,:)
      integer :: idx, nxc, nyc, rc
      integer,      allocatable :: itmp2d(:,:)
      real(kind=8), allocatable :: rtmp2d(:,:)
      character(128) :: msg=" "

      if ( present(olap) ) then
        local_olap = olap
      else
        local_olap = 0
      endif

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

      !--- Open the grid description file for reading
      ncstat = nf90_open(file_name, NF90_NOWRITE, ncfid)
      msg = "read_scrip_grid_file: nf90_open("//trim(file_name)//",..."
      call nc_error_handler(ncfid, ncstat, 1, msg=msg)

      !--- Verify that the source grid rank is 2
      ncstat = nf90_inq_dimid(ncfid, 'grid_rank', ncvid)
      msg = "read_scrip_grid_file: nf90_inq_dimid(...,grid_rank"
      call nc_error_handler(ncfid, ncstat, 10, msg=msg)
      ncstat = nf90_inquire_dimension(ncfid, ncvid, len=grid_rank)
      msg = "read_scrip_grid_file: nf90_inquire_dimension(..."
      call nc_error_handler(ncfid, ncstat, 11, msg=msg)
      if ( grid_rank /= 2 ) then
        write(*,*)'read_scrip_grid_file: Invalid grid rank = ',grid_rank
        call xit("read_scrip_grid_file",-1)
      endif

      !--- Verify that there are exactly 4 grid corners
      ncstat = nf90_inq_dimid(ncfid, 'grid_corners', ncvid)
      call nc_error_handler(ncfid, ncstat, 10)
      ncstat = nf90_inquire_dimension(ncfid, ncvid, len=grid_corners)
      call nc_error_handler(ncfid, ncstat, 11)
      if ( grid_corners /= 4 ) then
        write(*,*)'read_scrip_grid_file: Invalid grid corners = ',grid_corners
        call xit("read_scrip_grid_file",-2)
      endif

      !--- Verify that the number of grid edges is 2
      ncstat = nf90_inq_dimid(ncfid, 'edges', ncvid)
      call nc_error_handler(ncfid, ncstat, 10)
      ncstat = nf90_inquire_dimension(ncfid, ncvid, len=edges)
      call nc_error_handler(ncfid, ncstat, 11)
      if ( edges /= 2 ) then
        write(*,*)'read_scrip_grid_file: Invalid number of edges = ',edges
        call xit("read_scrip_grid_file",-1)
      endif

      if ( present(title) ) then
        !--- Read the grid "title" attribute
        local_title = ' '
        ncstat = nf90_get_att(ncfid,NF90_GLOBAL,'title',local_title)
        call nc_error_handler(ncfid, ncstat, 2)
        title = trim(local_title)
      endif

      !--- Always read the shape of the grid
      ncstat = nf90_inq_varid(ncfid, 'grid_dims', ncvid)
      call nc_error_handler(ncfid, ncstat, 18)
      ncstat = nf90_get_var(ncfid, ncvid, local_grid_dims)
      call nc_error_handler(ncfid, ncstat, 19)
      nx = local_grid_dims(1)
      ny = local_grid_dims(2)

      if ( present(grid_dims) ) then
        grid_dims(1:2) = local_grid_dims(1:2)
      endif

      !--- Always read the number of overlap grid points
      allocate( local_grid_edges(edges,grid_rank) )
      ncstat = nf90_inq_varid(ncfid, 'grid_edges', ncvid)
      call nc_error_handler(ncfid, ncstat, 18)
      ncstat = nf90_get_var(ncfid, ncvid, local_grid_edges)
      call nc_error_handler(ncfid, ncstat, 19)

      if ( present(grid_edges) ) then
        if ( size(grid_edges) < size(local_grid_edges) ) then
          write(6,*)"read_scrip_grid_file: User supplied grid_edges array is too small."
          write(6,*)"    size(grid_edges) = ",size(grid_edges)
          write(6,*)"       required size = ",size(local_grid_edges)
          call xit("read_scrip_grid_file",-1)
        endif
        grid_edges = local_grid_edges
      endif

      if ( .not. present(olap) .and. .not. keep_overlap ) then
        !--- Reset local_olap unless the user has supplied a value for olap
        !--- This implies that overlaping longitude columns will always be removed
        !--- when present unless the user supplies either olap=0 or overlap=.true.
        local_olap = local_grid_edges(2,1)
      endif
      deallocate( local_grid_edges )

      if ( local_olap < 0 .or. local_olap > nx-1 ) then
        !--- local_olap is the number of overlapping longitudes (coord dimension 1)
        !--- that will be removed from all fields returned by this routine
        write(6,*)"read_scrip_grid_file: olap = ",local_olap," is out of range."
        call xit("read_scrip_grid_file",-2)
      endif

      !--- Always read the grid size
      ncstat = nf90_inq_dimid(ncfid, 'grid_size', ncvid)
      call nc_error_handler(ncfid, ncstat, 6)
      ncstat = nf90_inquire_dimension(ncfid, ncvid, len=local_grid_size)
      call nc_error_handler(ncfid, ncstat, 7)

      if ( present(grid_size) ) then
        grid_size = local_grid_size
      endif

      if ( local_grid_size /= nx*ny ) then
        write(6,'(a)')"read_scrip_grid_file: Shape does not match grid size."
        write(6,'(a,3i8)')"       nx, ny, grid_size: ",nx,ny,local_grid_size
        call xit("read_scrip_grid_file",-3)
      endif

      if ( present(imask) ) then
        ncstat = nf90_inq_varid(ncfid, 'grid_imask', ncvid)
        call nc_error_handler(ncfid, ncstat, 30)
        if ( allocated(grid_imask) ) deallocate( grid_imask )
        allocate( grid_imask(local_grid_size) )
        ncstat = nf90_get_var(ncfid, ncvid, grid_imask)
        call nc_error_handler(ncfid, ncstat, 31)

        if ( associated(imask) ) nullify( imask )
        allocate( imask(nx,ny) )
        imask = reshape( grid_imask, (/nx,ny/) )
        deallocate( grid_imask )

        if ( local_olap > 0 .and. local_olap < nx ) then
          !--- Remove overlap longitudes
          !--- NOTE: No checks to verify these extra longitudes exist
          allocate( itmp2d(nx-local_olap,ny) )
          itmp2d = imask(1:nx-local_olap,:)
          deallocate( imask )
          allocate( imask(nx-local_olap,ny) )
          imask = itmp2d
          deallocate( itmp2d )
        endif

        if ( present(use_mask) ) then
          !--- use_mask = False means imask is set to 1 everywhere
          !--- use_mask = True  means imask is what was read from the file
          if ( .not. use_mask ) then
            imask = 1
            write(6,*)"read_scrip_grid_file: Setting imask = 1 over entire grid."
            call flush(6)
          endif
        endif
      endif

      if ( present(mid_lon) .or. present(mid_lon_1d) ) then
        ncstat = nf90_inq_varid(ncfid, 'grid_center_lon', ncvid)
        call nc_error_handler(ncfid, ncstat, 32)
        if ( allocated(grid_center_lon) ) deallocate( grid_center_lon )
        allocate( grid_center_lon(local_grid_size) )
        ncstat = nf90_get_var(ncfid, ncvid, grid_center_lon)
        call nc_error_handler(ncfid, ncstat, 33)

        !--- Assign grid_center_lon to a temporary 2d array
        if ( allocated(rtmp2d) ) deallocate(rtmp2d)
        allocate( rtmp2d(nx,ny) )
        rtmp2d = reshape( grid_center_lon, (/nx,ny/) )

        if ( present(mid_lon) ) then
          if ( associated(mid_lon) ) nullify( mid_lon )
          !--- Remove overlap longitudes when local_olap > 0
          !--- NOTE: No checks to verify these extra longitudes exist
          allocate( mid_lon(nx-local_olap,ny) )
          mid_lon = rtmp2d(1:nx-local_olap,1:ny)
        endif

        if ( present(mid_lon_1d) ) then
          if ( associated(mid_lon_1d) ) nullify( mid_lon_1d )
          !--- Remove overlap longitudes when local_olap > 0
          !--- NOTE: No checks to verify these extra longitudes exist
          allocate( mid_lon_1d(nx-local_olap) )
          mid_lon_1d = rtmp2d(1:nx-local_olap,1)
        endif
        deallocate( grid_center_lon )
        deallocate( rtmp2d )
      endif

      if ( present(mid_lat) .or. present(mid_lat_1d)  ) then
        ncstat = nf90_inq_varid(ncfid, 'grid_center_lat', ncvid)
        call nc_error_handler(ncfid, ncstat, 34)
        if ( allocated(grid_center_lat) ) deallocate( grid_center_lat )
        allocate( grid_center_lat(local_grid_size) )
        ncstat = nf90_get_var(ncfid, ncvid, grid_center_lat)
        call nc_error_handler(ncfid, ncstat, 35)

        !--- Assign grid_center_lat to a temporary 2d array
        if ( allocated(rtmp2d) ) deallocate(rtmp2d)
        allocate( rtmp2d(nx,ny) )
        rtmp2d = reshape( grid_center_lat, (/nx,ny/) )

        if ( present(mid_lat) ) then
          if ( associated(mid_lat) ) nullify( mid_lat )
          !--- Remove overlap longitudes when local_olap > 0
          !--- NOTE: No checks to verify these extra longitudes exist
          allocate( mid_lat(nx-local_olap,ny) )
          mid_lat = rtmp2d(1:nx-local_olap,1:ny)
        endif

        if ( present(mid_lat_1d) ) then
          if ( associated(mid_lat_1d) ) nullify( mid_lat_1d )
          !--- Remove overlap longitudes when local_olap > 0
          !--- NOTE: No checks to verify these extra longitudes exist
          allocate( mid_lat_1d(ny) )
          mid_lat_1d = rtmp2d(1,1:ny)
        endif
        deallocate( grid_center_lat )
        deallocate( rtmp2d )
      endif

      if ( (present(corn_lon) .and. present(corn_lat)) .or. &
           (present(corn_lon_1d) .and. present(corn_lat_1d)) ) then
        !--- Read corner lon from the file
        ncstat = nf90_inq_varid(ncfid, 'grid_corner_lon', ncvid)
        call nc_error_handler(ncfid, ncstat, 36)
        if ( allocated(grid_corner_lon) ) deallocate( grid_corner_lon )
        allocate( grid_corner_lon(grid_corners, local_grid_size) )
        ncstat = nf90_get_var(ncfid, ncvid, grid_corner_lon)
        call nc_error_handler(ncfid, ncstat, 37)

        !--- Read corner lat from the file
        ncstat = nf90_inq_varid(ncfid, 'grid_corner_lat', ncvid)
        call nc_error_handler(ncfid, ncstat, 38)
        if ( allocated(grid_corner_lat) ) deallocate( grid_corner_lat )
        allocate( grid_corner_lat(grid_corners, local_grid_size) )
        ncstat = nf90_get_var(ncfid, ncvid, grid_corner_lat)
        call nc_error_handler(ncfid, ncstat, 39)

        !--- Put grid cell corners into 2D arrays dimensioned (nx+1,ny+1)
        if ( allocated(corn_lon_tmp) ) deallocate(corn_lon_tmp)
        allocate( corn_lon_tmp(nx+1,ny+1) )
        if ( allocated(corn_lat_tmp) ) deallocate(corn_lat_tmp)
        allocate( corn_lat_tmp(nx+1,ny+1) )
        call convert_corners(nx, ny, grid_corner_lon, grid_corner_lat, &
                                     corn_lon_tmp,    corn_lat_tmp)

        if ( present(corn_lon) ) then
          if ( associated(corn_lon) ) nullify( corn_lon )
          !--- Remove overlap longitudes when local_olap > 0
          !--- NOTE: No checks to verify these extra longitudes exist
          allocate( corn_lon(nx+1-local_olap,ny+1) )
          corn_lon = corn_lon_tmp(1:nx+1-local_olap,1:ny+1)
        endif

        if ( present(corn_lat) ) then
          if ( associated(corn_lat) ) nullify( corn_lat )
          !--- Remove overlap longitudes when local_olap > 0
          !--- NOTE: No checks to verify these extra longitudes exist
          allocate( corn_lat(nx+1-local_olap,ny+1) )
          corn_lat = corn_lat_tmp(1:nx+1-local_olap,1:ny+1)
        endif

        if ( present(corn_lon_1d) ) then
          if ( associated(corn_lon_1d) ) nullify( corn_lon_1d )
          !--- Remove overlap longitudes when local_olap > 0
          !--- NOTE: No checks to verify these extra longitudes exist
          allocate( corn_lon_1d(nx+1-local_olap) )
          corn_lon_1d = corn_lon_tmp(1:nx+1-local_olap,1)
        endif

        if ( present(corn_lat_1d) ) then
          if ( associated(corn_lat_1d) ) nullify( corn_lat_1d )
          !--- Remove overlap longitudes when local_olap > 0
          !--- NOTE: No checks to verify these extra longitudes exist
          allocate( corn_lat_1d(ny+1) )
          corn_lat_1d = corn_lat_tmp(1,1:ny+1)
        endif

        deallocate( corn_lon_tmp, corn_lat_tmp )
        deallocate( grid_corner_lon, grid_corner_lat )

      else if ( present(corn_lon) .and. .not. present(corn_lat) ) then
        write(6,*)"read_scrip_grid_file: corn_lon is present without corn_lat"
        call xit("read_scrip_grid_file", -1)

      else if ( present(corn_lat) .and. .not. present(corn_lon) ) then
        write(6,*)"read_scrip_grid_file: corn_lat is present without corn_lon"
        call xit("read_scrip_grid_file", -1)

      else if ( present(corn_lon_1d) .and. .not. present(corn_lat_1d) ) then
        write(6,*)"read_scrip_grid_file: corn_lon_1d is present without corn_lat_1d"
        call xit("read_scrip_grid_file", -1)

      else if ( present(corn_lat_1d) .and. .not. present(corn_lon_1d) ) then
        write(6,*)"read_scrip_grid_file: corn_lat_1d is present without corn_lon_1d"
        call xit("read_scrip_grid_file", -1)

      endif

      if ( present(grid_dims) ) then
        !--- If local_olap > 0 then reset grid_dims(1) to be consistent
        !--- with the shape of the 2D arrays that are returned
        grid_dims(1) = grid_dims(1) - local_olap
      endif

      !--- Close the grid description file
      ncstat = nf90_close(ncfid)

    end subroutine read_scrip_grid_file

    !***************************************************************************
    !--- Convert 2D corner coordinate arrays dimensioned (4, nx*ny)
    !--- into a 2D arrays dimensioned (nx+1, ny+1)
    !--- This assumes that there are exactly 4 corners per grid cell and all
    !--- corners are ordered the same (counterclockwise) around each cell center
    !--- This also assumes that the cells are layed out in column major order
    !--- along the second dimension of input arrays corner_lon and corner_lat
    !***************************************************************************
    subroutine convert_corners(nx, ny, corner_lon,   corner_lat, &
                                       corn_lon_out, corn_lat_out)
      integer(kind=4) :: nx, ny
      !--- corner_lon and corner_lat must be dimensioned (4,nx*ny)
      real(kind=8)    :: corner_lon(:,:), corner_lat(:,:)
      !--- corn_lon_out and corn_lat_out must be dimensioned (nx+1,ny+1)
      real(kind=8)    :: corn_lon_out(:,:), corn_lat_out(:,:)

      !--- Local
      integer(kind=4) :: indx, ix, iy

      if ( size(corner_lon,dim=1) /= 4     .or. &
           size(corner_lon,dim=2) /= nx*ny ) then
        write(6,*)"convert_corners: corner_lon shape mismatch."
        write(6,*)"    shape(corner_lon) = ",shape(corner_lon)
        write(6,*)"              expecting   4 ",nx*ny
        call xit("convert_corners", -1)
      endif

      if ( size(corner_lat,dim=1) /= 4     .or. &
           size(corner_lat,dim=2) /= nx*ny ) then
        write(6,*)"convert_corners: corner_lat shape mismatch."
        write(6,*)"    shape(corner_lat) = ",shape(corner_lat)
        write(6,*)"              expecting   4 ",nx*ny
        call xit("convert_corners", -1)
      endif

      if ( size(corn_lon_out,dim=1) /= nx+1 .or. &
           size(corn_lon_out,dim=2) /= ny+1 ) then
        write(6,*)"convert_corners: corn_lon_out shape mismatch."
        write(6,*)"    shape(corn_lon_out) = ",shape(corn_lon_out)
        write(6,*)"              expecting   ",nx+1,ny+1
        call xit("convert_corners", -1)
      endif

      if ( size(corn_lat_out,dim=1) /= nx+1 .or. &
           size(corn_lat_out,dim=2) /= ny+1 ) then
        write(6,*)"convert_corners: corn_lat_out shape mismatch."
        write(6,*)"    shape(corn_lat_out) = ",shape(corn_lat_out)
        write(6,*)"              expecting   ",nx+1,ny+1
        call xit("convert_corners", -1)
      endif

      !--- Example 2x3 grid in column major order with corners numbered
      !--- Note that rows have index increasing from bottom to top
      !---  |-------|-------|-------|
      !---  |4     3|4     3|4     3|
      !---  |  2,1  |  2,2  |  2,3  |
      !---  |1     2|1     2|1     2|
      !---  |-------|-------|-------|
      !---  |4     3|4     3|4     3|
      !---  |  1,1  |  1,2  |  1,3  |
      !---  |1     2|1     2|1     2|
      !---  |-------|-------|-------|

      do indx=1,nx*ny
        !--- Determine indicies on the 2D grid assuming column major order
        ix = 1 + mod(indx-1_4,nx)
        iy = (indx-1)/nx + 1

        !--- Assign bottom left corner values to all cells as default
        corn_lon_out(ix,iy) = corner_lon(1,indx)
        corn_lat_out(ix,iy) = corner_lat(1,indx)

        if ( ix == nx ) then
          !--- This is the last grid cell in this row
          !--- Overwrite these cells with bottom right corner values
          corn_lon_out(ix+1,iy) = corner_lon(2,indx)
          corn_lat_out(ix+1,iy) = corner_lat(2,indx)
        endif

        if ( iy == ny ) then
          !--- This is the last grid cell in this column
          !--- Overwrite these cells with top left corner values
          corn_lon_out(ix,iy+1) = corner_lon(4,indx)
          corn_lat_out(ix,iy+1) = corner_lat(4,indx)
        endif
      enddo

      !--- This is the top right corner cell in the grid
      !--- Overwrite this cell with the top right corner value
      corn_lon_out(nx+1,ny+1) = corner_lon(3,nx*ny)
      corn_lat_out(nx+1,ny+1) = corner_lat(3,nx*ny)

    end subroutine convert_corners

    !***************************************************************************
    !--- Read AGCM grid cell area from a file
    !***************************************************************************
    subroutine read_agcm_grid_area_file(file_name, area)

      implicit none

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

      !--- The grid cell area on the AGCM grid (m2)
      real(kind=8), intent(out), pointer :: area(:,:)

      !--- Local
      logical :: exists
      integer(kind=4) :: fid

      if ( len_trim(file_name) < 1 ) then
        write(6,*)"read_agcm_grid_area_file: file_name is blank"
        call flush(6)
        call xit("read_agcm_grid_area_file",-1)
      endif

      inquire(file=trim(file_name), exist=exists)
      if ( .not. exists ) then
        write(6,*)"read_agcm_grid_area_file: File does not exist ",trim(file_name)
        call flush(6)
        call xit("read_agcm_grid_area_file",-2)
      endif

      !--- Open the grid area file for reading
      call ncdf_open(fid, trim(file_name), mode="read")

      !--- Read area from this file
      call ncdf_read_var(fid, "areacella", area)

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

    end subroutine read_agcm_grid_area_file

    !***************************************************************************
    !--- Read specific data from a NEMO "mesh_mask" file
    !***************************************************************************
    subroutine read_NEMO_mesh_mask_file(file_name, grid, overlap, &
                                        mid_lon, mid_lat, corn_lon, corn_lat, &
                                        mask, imask, area, &
                                        title, grid_size, grid_dims, grid_edges)

      implicit none

      !--- The name of the NEMO mesh_mask file
      character(*), intent(in)  :: file_name

      !--- The NEMO grid type to extract (T,U,V,F) default is T grid
      character(*), intent(in), optional :: grid

      !--- If overlap = F then do not include overlaping grid cells in any
      !--- dimension in the output grid description
      logical, intent(in), optional :: overlap

      !--- lon/lat of grid cell centers
      real(kind=8), intent(out), pointer, optional ::  mid_lon(:,:)
      real(kind=8), intent(out), pointer, optional ::  mid_lat(:,:)

      !--- lon/lat of grid cell corners
      !--- It is assumed that there are 4 corners per cell
      real(kind=8), intent(out), pointer, optional :: corn_lon(:,:)
      real(kind=8), intent(out), pointer, optional :: corn_lat(:,:)

      !--- Land mask at level 1 of the T grid (real or integer)
      real(kind=8), intent(out), pointer, optional ::  mask(:,:)
      integer,      intent(out), pointer, optional :: imask(:,:)

      !--- The grid cell area (m2)
      real(kind=8), intent(out), pointer, optional :: area(:,:)

      !--- The grid title
      character(*), intent(out), optional :: title

      !--- The total number of points in the grid
      integer, intent(out), optional :: grid_size

      !--- The shape of the grid ie (nx,ny)
      integer, intent(out), optional :: grid_dims(2)

      !--- The number of overlap grids (lower and upper) in each dimension
      integer, intent(out), optional :: grid_edges(2,2)

      !--- Local
      logical :: exists, resize
      integer(kind=4) :: ncstat, fid, varid
      integer(kind=4) :: ndims, dimids(7), vstart(7), vcount(7)
      character(256) :: strng
      character(128) :: msg=" "
      character(32)  :: lon_name="glamt", lat_name="gphit"
      character(32)  :: lon_corn_name="glamf", lat_corn_name="gphif"
      character(32)  :: e1_name="e1t", e2_name="e2t", mask_name="tmask"
      character(1)   :: gtype="T"
      real(kind=8)   :: gsize
      real(kind=8), pointer :: awrk2d(:,:), bwrk2d(:,:)
      integer,      pointer :: iwrk2d(:,:)
      real(kind=8), pointer, dimension(:,:)   :: clon=>null(), clat=>null()
      real(kind=8), pointer, dimension(:,:)   :: e1=>null(), e2=>null()
      real(kind=8), pointer, dimension(:,:,:) :: xmask=>null()
      integer     , pointer, dimension(:,:,:) :: ximask=>null()
      integer :: nx, ny, idx, ic, jc
      integer :: offset(2)
      character(64),   pointer :: dim_names(:)
      integer,         pointer :: dim_sizes(:)

      if ( len_trim(file_name) < 1 ) then
        write(6,*)"read_NEMO_mesh_mask_file: file_name is blank"
        call flush(6)
        call xit("read_NEMO_mesh_mask_file",-1)
      endif

      inquire(file=trim(file_name), exist=exists)
      if ( .not. exists ) then
        write(6,*)"read_NEMO_mesh_mask_file: File does not exist ",trim(file_name)
        call flush(6)
        call xit("read_NEMO_mesh_mask_file",-2)
      endif

      if ( present(overlap) ) then
        !--- If overlap = T then keep the overlap longitude in output arrays
        !--- If overlap = F then remove the overlap longitude before return
        resize = .not. overlap
      else
        !--- Remove the overlap longitude by default
        resize = .true.
      endif

      if ( present(grid) ) then
        gtype = grid(1:1)
        !--- Define names of variables found in the mesh mask file that correspond to
        !--- grid cell centers, corners and scale factors
        !--- offset is what needs to be added to (i,j) to get to the correct elements
        !--- in the corner arrays (as read from the file) that correspond to the
        !--- northeast corner of each cell
        select case(gtype)
          case ("t", "T")
            gtype = 'T'
            lon_name = "glamt"
            lat_name = "gphit"
            lon_corn_name = "glamf"
            lat_corn_name = "gphif"
            mask_name = "tmask"
            e1_name = "e1t"
            e2_name = "e2t"
            offset = (/0,0/)
          case ("u", "U")
            gtype = 'U'
            lon_name = "glamu"
            lat_name = "gphiu"
            lon_corn_name = "glamv"
            lat_corn_name = "gphiv"
            mask_name = "umask"
            e1_name = "e1u"
            e2_name = "e2u"
            offset = (/1,0/)
          case ("v", "V")
            gtype = 'V'
            lon_name = "glamv"
            lat_name = "gphiv"
            lon_corn_name = "glamu"
            lat_corn_name = "gphiu"
            mask_name = "vmask"
            e1_name = "e1v"
            e2_name = "e2v"
            offset = (/0,1/)
          case ("f", "F")
            !--- The values for lon_corn_name, lat_corn_name and the values for the offset
            !--- need to be verified for the F grid
            write(6,*)"read_NEMO_mesh_mask_file: Grid type F not currently supported"
            call flush(6)
            call xit("read_NEMO_mesh_mask_file",-9)
            gtype = 'F'
            lon_name = "glamf"
            lat_name = "gphif"
            lon_corn_name = "glamt"
            lat_corn_name = "gphit"
            mask_name = "fmask"
            e1_name = "e1f"
            e2_name = "e2f"
            offset = (/-1,0/)
          case default
            write(6,*)"read_NEMO_mesh_mask_file: Invalid grid type ",gtype
            call flush(6)
            call xit("read_NEMO_mesh_mask_file",-3)
        end select
      else
        !--- Default is the T grid
        gtype = 'T'
        lon_name = "glamt"
        lat_name = "gphit"
        lon_corn_name = "glamf"
        lat_corn_name = "gphif"
        mask_name = "tmask"
        e1_name = "e1t"
        e2_name = "e2t"
        offset = (/0,0/)
      endif

      !--- Open the mesh_mask file for reading
      call ncdf_open(fid, trim(file_name), mode="read")

      !--- Read dimension sizes from the file
      call ncdf_inquire_dims(fid, dim_names=dim_names, dim_sizes=dim_sizes)
      nx = -1
      ny = -1
      do idx=1,size(dim_names)
        if ( trim(adjustl(dim_names(idx))) .eq. "x") nx = dim_sizes(idx)
        if ( trim(adjustl(dim_names(idx))) .eq. "y") ny = dim_sizes(idx)
      enddo
      if ( nx < 0 .or. ny < 0 ) then
        write(6,*)"read_NEMO_mesh_mask_file: Unable to determine nx/ny from ",trim(file_name)
        call flush(6)
        call xit("read_NEMO_mesh_mask_file",-4)
      endif

      write(6,*)"read_NEMO_mesh_mask_file: ",trim(file_name),"   nx=",nx,"   ny=",ny
      call flush(6)

      if ( present(title) ) then
        title = " "
        if ( nx > 999 .or. ny > 999 ) then
          write(title,'(3a,i4,a,i4)')"ORCA ",gtype," grid -- ",nx,"x",ny
        else
          write(title,'(3a,i3.3,a,i3.3)')"ORCA ",gtype," grid -- ",nx,"x",ny
        endif
        write(6,*)"read_NEMO_mesh_mask_file: title = ",trim(title)
        call flush(6)
      endif

      if ( present(grid_edges) ) then
        grid_edges = 0
        !--- Assume 2 overlap logitudes
        grid_edges(2,1) = 2
        write(6,*)"read_NEMO_mesh_mask_file: grid_edges = ",grid_edges
        call flush(6)
      endif

      if ( present(mid_lon) ) then
        !--- The file should contain
        !---     float glamt(t, y, x)
        !---     float glamu(t, y, x)
        !---     float glamv(t, y, x)
        !---     float glamf(t, y, x)
        !--- Read grid cell center longitudes from the file
        !--- This call will destroy then reallocate mid_lon
        call ncdf_read_var(fid, trim(lon_name), mid_lon, nt=1)

        if ( resize ) then
          !--- Remove 2 overlap longitudes
          nx = size(mid_lon,dim=1) - 2
          ny = size(mid_lon,dim=2)
!          call resize_a2d(mid_lon,nx,ny)
          call resize_nemo_field(mid_lon)
        endif
      endif

      if ( present(mid_lat) ) then
        !--- The file should contain
        !---     float gphit(t, y, x)
        !---     float gphiu(t, y, x)
        !---     float gphiv(t, y, x)
        !---     float gphif(t, y, x)
        !--- Read grid cell center latitudes from the file
        !--- This call will destroy then reallocate mid_lat
        call ncdf_read_var(fid, trim(lat_name), mid_lat, nt=1)

        if ( resize ) then
          !--- Remove 2 overlap longitudes
          nx = size(mid_lat,dim=1) - 2
          ny = size(mid_lat,dim=2)
!          call resize_a2d(mid_lat,nx,ny)
          call resize_nemo_field(mid_lat)
        endif
      endif

      if ( present(corn_lon) ) then
        !--- The file should contain
        !---     float glamt(t, y, x)
        !---     float glamu(t, y, x)
        !---     float glamv(t, y, x)
        !---     float glamf(t, y, x)
        !--- Read grid cell corner longitudes from the file
        !--- This call will destroy then reallocate clon
        call ncdf_read_var(fid, trim(lon_corn_name), clon, nt=1)

        !--- Reformat these corner lon values to add missing boundary points
        !--- corn_lon is returned with 1 more row and 1 more column than clon
        call mm_corner_lon(clon, offset, corn_lon)

        if ( resize ) then
          !--- Remove 2 overlap longitudes
          nx = size(corn_lon,dim=1) - 2
          ny = size(corn_lon,dim=2)
!          call resize_a2d(corn_lon,nx,ny)
          call resize_nemo_field(corn_lon,corner=.true.)
        endif
      endif

      if ( present(corn_lat) ) then
        !--- The file should contain
        !---     float gphit(t, y, x)
        !---     float gphiu(t, y, x)
        !---     float gphiv(t, y, x)
        !---     float gphif(t, y, x)
        !--- Read grid cell corner latitudes from the file
        !--- This call will destroy then reallocate clat
        call ncdf_read_var(fid, trim(lat_corn_name), clat, nt=1)

        !--- Reformat these corner lat values to add missing boundary points
        !--- corn_lat is returned with 1 more row and 1 more column than clat
        call mm_corner_lat(clat, offset, corn_lat)

        if ( resize ) then
          !--- Remove 2 overlap longitudes
          nx = size(corn_lat,dim=1) - 2
          ny = size(corn_lat,dim=2)
!          call resize_a2d(corn_lat,nx,ny)
          call resize_nemo_field(corn_lat,corner=.true.)
        endif
      endif

      if ( present(mask) ) then
!--- TODO --- This should return tmask_i
write(6,*)"read_NEMO_mesh_mask_file: Optional mask argument not yet supported"
call xit("read_NEMO_mesh_mask_file",-98)
        !--- The file should contain
        !---     byte tmask(t, z, y, x)
        !---     byte umask(t, z, y, x)
        !---     byte vmask(t, z, y, x)
        !---     byte fmask(t, z, y, x)
        !--- Read T grid land mask from the file
        call ncdf_read_var(fid, trim(mask_name), xmask, nt=1)

        !--- Assign mask with surface (level 1) mask values from
        !--- the 3D mask that was read from the file
        if ( associated(mask) ) deallocate(mask)
        allocate( mask(size(xmask,dim=1),size(xmask,dim=2)) )
        mask = xmask(:,:,1)
        deallocate(xmask)

        if ( resize ) then
          !--- Remove 2 overlap longitudes
          nx = size(mask,dim=1) - 2
          ny = size(mask,dim=2)
!          call resize_a2d(mask,nx,ny)
          call resize_nemo_field(mask)
        endif
      endif

      if ( present(imask) ) then
!--- TODO --- This should return tmask_i
write(6,*)"read_NEMO_mesh_mask_file: Optional imask argument not yet supported"
call xit("read_NEMO_mesh_mask_file",-99)
        !--- The file should contain
        !---     byte tmask(t, z, y, x)
        !---     byte umask(t, z, y, x)
        !---     byte vmask(t, z, y, x)
        !---     byte fmask(t, z, y, x)
        !--- Read T grid land mask from the file
        call ncdf_read_var(fid, trim(mask_name), ximask, nt=1)

        !--- Assign mask with surface (level 1) mask values from
        !--- the 3D mask that was read from the file
        if ( associated(imask) ) deallocate(imask)
        allocate( imask(size(ximask,dim=1),size(ximask,dim=2)) )
        imask = ximask(:,:,1)
        deallocate(ximask)

        if ( resize ) then
          !--- Remove 2 overlap longitudes
          nx = size(imask,dim=1) - 2
          ny = size(imask,dim=2)
!          call resize_a2d(imask,nx,ny)
          call resize_nemo_field(imask)
        endif
      endif

      if ( present(area) ) then
        !--- The file should contain
        !---     double e1t(t, y, x)
        !---     double e1u(t, y, x)
        !---     double e1v(t, y, x)
        !---     double e1f(t, y, x)
        !---
        !---     double e2t(t, y, x)
        !---     double e2u(t, y, x)
        !---     double e2v(t, y, x)
        !---     double e2f(t, y, x)
        !--- Read scale lengths from the file

        call ncdf_read_var(fid, trim(e1_name), e1, nt=1)
        call ncdf_read_var(fid, trim(e2_name), e2, nt=1)

        !--- Grid cell area is the product of these two scale lengths
        if ( associated(area) ) deallocate(area)
        allocate( area(size(e1,dim=1),size(e1,dim=2)) )
        area = e1 * e2
        deallocate( e1, e2 )

        if ( resize ) then
          !--- Remove 2 overlap longitudes
          nx = size(area,dim=1) - 2
          ny = size(area,dim=2)
!          call resize_a2d(area,nx,ny)
          call resize_nemo_field(area)
        endif

      endif

      !--- Define grid_dims and grid_size here so that any changes to nx,ny due to
      !--- removing overlap longitudes will be reflected in these output values
      if ( present(grid_dims) ) then
        grid_dims(1) = nx
        grid_dims(2) = ny
      endif

      if ( present(grid_size) ) then
        grid_size = nx*ny
      endif

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

    end subroutine read_NEMO_mesh_mask_file

    subroutine mm_corner_lon(clon, offset, corn_lon)
      !--- clon contains corner longitudes as read from the mesh_mask file
      real(kind=8), intent(in),  pointer :: clon(:,:)

      !--- offset is what needs to be added to (i,j) to get to the correct elements
      !--- in the clon arrays (read from the mesh_mask file) that correspond to the
      !--- northeast corner of each cell
      integer, intent(in) :: offset(2)

      !--- corn_lon is output and will contain the corner longitudes
      !--- suitable for use in ESMF grid definitions
      real(kind=8), intent(out), pointer :: corn_lon(:,:)

      !--- Local
      real(kind=8), pointer, dimension(:,:,:) :: corner_lon=>null()
      real(kind=8), parameter :: circle = 360.0_8 ! assumes corn_lon is in degrees
      integer :: nx, ny, idx, n, i, j, ic, jc, imid, isame

      if ( .not. associated(clon) ) then
        write(6,*)"mm_corner_lon: Input clon array is not allocated"
        call flush(6)
        call xit("mm_corner_lon",-1)
      endif

      !--- Get the dimension size from the input clon array
      nx = size(clon,dim=1)
      ny = size(clon,dim=2)

      !--- Always destroy then recreate the corn_lon output array
      if ( associated(corn_lon) ) deallocate(corn_lon)
      allocate ( corn_lon(nx+1,ny+1) )
      corn_lon = 0.0_8
  
      !-----------------------------------------------------------------------
      ! Corners must be ordered as follows:
      !     4 3
      !     1 2
      !
      ! Create a pair of arrays for the corners where corn_lon(1,1) corresponds
      ! to the south west corner of a box centered at (glam(1,1),gphi(1,1))

      !--- Assign the internal points
      do j = 2,ny
        do i = 2,nx
          ic = i + offset(1) - 1
          jc = j + offset(2) - 1
          corn_lon(i,j) = clon(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
      !---     clon(1,:) == clon(nx-1,:)   clat(1,:) == clat(nx-1,:)
      !---     clon(2,:) == clon(nx,:)     clat(2,:) == clat(nx,:)
      do i=2,nx
        ic = i + offset(1) - 1
        !--- First (south/bottom) row
        corn_lon(i,1) = clon(ic,1)

        !--- Last (north/top) row
        corn_lon(i,ny+1) = clon(ic,ny)
      enddo

      do j=2,ny
        jc = j + offset(2) - 1
        !--- First (west/left) column
        corn_lon(1,j) = clon(nx-2,jc)

        !--- Last (east/right) column
        corn_lon(nx+1,j) = clon(2,jc)
      enddo

      !--- Then the corner points
      ic = nx-2
      corn_lon(1,1) = clon(ic,1)

      ic = 2
      corn_lon(nx+1,1) = clon(ic,1)

      imid = (nx-1)/2 + 1
      ic = 2*imid - nx + 1
      corn_lon(1,ny+1) = clon(ic,ny)

      ic = 2*imid
      corn_lon(nx+1,ny+1) = clon(ic,ny)

      !--- Zero out small values
      where (abs(corn_lon)<1e-20) corn_lon=0.0_8

    end subroutine mm_corner_lon

    subroutine mm_corner_lat(clat, offset, corn_lat)
      !--- clat contains corner latitudes as read from the mesh_mask file
      real(kind=8), intent(in),  pointer :: clat(:,:)

      !--- offset is what needs to be added to (i,j) to get to the correct elements
      !--- in the clat arrays (read from the mesh_mask file) that correspond to the
      !--- northeast corner of each cell
      integer, intent(in) :: offset(2)

      !--- corn_lat is output and will contain the corner latitudes
      !--- suitable for use in ESMF grid definitions
      real(kind=8), intent(out), pointer :: corn_lat(:,:)

      !--- Local
      real(kind=8), pointer, dimension(:,:,:) :: corner_lat=>null()
      integer :: nx, ny, idx, n, i, j, ic, jc, imid, isame

      if ( .not. associated(clat) ) then
        write(6,*)"mm_corner_lat: Input clat array is not allocated"
        call flush(6)
        call xit("mm_corner_lat",-1)
      endif

      !--- Get the dimension size from the input clat array
      nx = size(clat,dim=1)
      ny = size(clat,dim=2)

      !--- Always destroy then recreate the corn_lat output array
      if ( associated(corn_lat) ) deallocate(corn_lat)
      allocate ( corn_lat(nx+1,ny+1) )
      corn_lat = 0.0_8
  
      !-----------------------------------------------------------------------
      ! Corners must be ordered as follows:
      !     4 3
      !     1 2
      !
      ! Create a pair of arrays for the corners where corn_lat(1,1) corresponds
      ! to the south west corner of a box centered at (glam(1,1),gphi(1,1))
 
      !--- Assign the internal points
      do j = 2,ny
        do i = 2,nx
          ic = i + offset(1) - 1
          jc = j + offset(2) - 1
          corn_lat(i,j) = clat(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
      !---     clon(1,:) == clon(nx-1,:)   clat(1,:) == clat(nx-1,:)
      !---     clon(2,:) == clon(nx,:)     clat(2,:) == clat(nx,:)
      do i=2,nx
        ic = i + offset(1) - 1
        !--- First (south/bottom) row
        corn_lat(i,1) = clat(ic,1) - (clat(ic,2)-clat(ic,1))

        !--- Last (north/top) row
        corn_lat(i,ny+1) = clat(ic,ny)
      enddo

      do j=2,ny
        jc = j + offset(2) - 1
        !--- First (west/left) column
        corn_lat(1,j) = clat(nx-2,jc)

        !--- Last (east/right) column
        corn_lat(nx+1,j) = clat(2,jc)
      enddo

      !--- Then the corner points
      ic = nx-2
      corn_lat(1,1) = clat(ic,1) - (clat(ic,2)-clat(ic,1))

      ic = 2
      corn_lat(nx+1,1) = clat(ic,1) - (clat(ic,2)-clat(ic,1))

      imid = (nx-1)/2 + 1
      ic = 2*imid - nx + 1
      corn_lat(1,ny+1) = clat(ic,ny)

      ic = 2*imid
      corn_lat(nx+1,ny+1) = clat(ic,ny)

      !--- Zero out small values
      where (abs(corn_lat)<1e-20) corn_lat=0.0_8

    end subroutine mm_corner_lat

    !***************************************************************************
    !--- Resize the 2D input array to the input shape (nx,ny)
    !---     1 <= nx <= size(afld,dim=1)
    !---     1 <= ny <= size(afld,dim=2)
    !--- All rows greater than nx and/or columns greater than ny will be removed
    !***************************************************************************
    !***************************************************************************
    !--- Resize a real*4 input array
    !***************************************************************************
    subroutine resize_a2d_r4(afld,nx,ny)
      !--- real*4 input array
      real(kind=4), pointer, intent(inout) :: afld(:,:)
      integer, intent(in) :: nx, ny

      !--- Local
      real(kind=4), pointer, dimension(:,:) :: wrk2d=>null()

      if ( .not. associated(afld) ) then
        write(6,*)"resize_a2d_r4: Input array is not allocated."
        call xit("resize_a2d_r4",-1)
      endif

      if ( nx < 1 .or. ny < 1 .or. &
           nx > size(afld,dim=1) .or. ny > size(afld,dim=2) ) then
        write(6,*)"resize_a2d_r4: Invalid new shape"
        write(6,*)"   current shape: ",shape(afld)
        write(6,*)"       new shape: ",nx,ny
        call xit("resize_a2d_r4",-2)
      endif

      allocate( wrk2d(nx,ny) )
      wrk2d(1:nx,1:ny) = afld(1:nx,1:ny)

      deallocate( afld )
      allocate( afld(nx,ny) )
      afld(1:nx,1:ny) = wrk2d(1:nx,1:ny)

      deallocate( wrk2d )

    end subroutine resize_a2d_r4

    !***************************************************************************
    !--- Resize a real*8 input array
    !***************************************************************************
    subroutine resize_a2d_r8(afld,nx,ny)
      !--- real*8 input array
      real(kind=8), pointer, intent(inout) :: afld(:,:)
      integer, intent(in) :: nx, ny

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

      if ( .not. associated(afld) ) then
        write(6,*)"resize_a2d_r8: Input array is not allocated."
        call xit("resize_a2d_r8",-1)
      endif

      if ( nx < 1 .or. ny < 1 .or. &
           nx > size(afld,dim=1) .or. ny > size(afld,dim=2) ) then
        write(6,*)"resize_a2d_r8: Invalid new shape"
        write(6,*)"   current shape: ",shape(afld)
        write(6,*)"       new shape: ",nx,ny
        call xit("resize_a2d_r8",-2)
      endif

      allocate( wrk2d(nx,ny) )
      wrk2d(1:nx,1:ny) = afld(1:nx,1:ny)

      deallocate( afld )
      allocate( afld(nx,ny) )
      afld(1:nx,1:ny) = wrk2d(1:nx,1:ny)

      deallocate( wrk2d )

    end subroutine resize_a2d_r8

    !***************************************************************************
    !--- Resize an integer*4 input array
    !***************************************************************************
    subroutine resize_a2d_i4(afld,nx,ny)
      !--- integer*4 input array
      integer(kind=4), pointer, intent(inout) :: afld(:,:)
      integer, intent(in) :: nx, ny

      !--- Local
      integer(kind=4), pointer, dimension(:,:) :: wrk2d=>null()

      if ( .not. associated(afld) ) then
        write(6,*)"resize_a2d_i4: Input array is not allocated."
        call xit("resize_a2d_i4",-1)
      endif

      if ( nx < 1 .or. ny < 1 .or. &
           nx > size(afld,dim=1) .or. ny > size(afld,dim=2) ) then
        write(6,*)"resize_a2d_i4: Invalid new shape"
        write(6,*)"   current shape: ",shape(afld)
        write(6,*)"       new shape: ",nx,ny
        call xit("resize_a2d_i4",-2)
      endif

      allocate( wrk2d(nx,ny) )
      wrk2d(1:nx,1:ny) = afld(1:nx,1:ny)

      deallocate( afld )
      allocate( afld(nx,ny) )
      afld(1:nx,1:ny) = wrk2d(1:nx,1:ny)

      deallocate( wrk2d )

    end subroutine resize_a2d_i4

    !***************************************************************************
    !--- Resize an integer*8 input array
    !***************************************************************************
    subroutine resize_a2d_i8(afld,nx,ny)
      !--- integer*8 input array
      integer(kind=8), pointer, intent(inout) :: afld(:,:)
      integer, intent(in) :: nx, ny

      !--- Local
      integer(kind=8), pointer, dimension(:,:) :: wrk2d=>null()

      if ( .not. associated(afld) ) then
        write(6,*)"resize_a2d_i8: Input array is not allocated."
        call xit("resize_a2d_i8",-1)
      endif

      if ( nx < 1 .or. ny < 1 .or. &
           nx > size(afld,dim=1) .or. ny > size(afld,dim=2) ) then
        write(6,*)"resize_a2d_i8: Invalid new shape"
        write(6,*)"   current shape: ",shape(afld)
        write(6,*)"       new shape: ",nx,ny
        call xit("resize_a2d_i8",-2)
      endif

      allocate( wrk2d(nx,ny) )
      wrk2d(1:nx,1:ny) = afld(1:nx,1:ny)

      deallocate( afld )
      allocate( afld(nx,ny) )
      afld(1:nx,1:ny) = wrk2d(1:nx,1:ny)

      deallocate( wrk2d )

    end subroutine resize_a2d_i8

    !***************************************************************************
    !--- Add meta data for a given variable to a netcdf file
    !--- The variable must already be defined and the netcdf file must be in
    !--- define mode when this routine is called
    !***************************************************************************
    subroutine nc_add_var_meta_data(var_name, fid, vid)
      character(*),    intent(in) :: var_name
      integer(kind=4), intent(in) :: fid
      integer(kind=4), intent(in) :: vid

      !--- Local
      integer(kind=4) :: ncstat
      character(128)  :: strng=" ", strng2=" "
      character(32)   :: sn_snd_thick_cldes

      strng2 = "nc_add_var_meta_data: var="//trim(var_name)

      !--- This is a temporary placeholder to indicate what the name/units for
      !--- OIceTck and OSnwTck should be. This meta data could change when CICE
      !--- is used for the ice model.
      sn_snd_thick_cldes = "weighted iwe and swe"

      !--- Define attributes for known variable names
      !--- If the name is not known then do nothing
      select case (trim(adjustl(var_name)))
        case ("GT_atm")
          !--- Ground temperature
          strng = "Ground Temperature"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "K"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("GTICE_atm")
          !--- Sea Ice Surface temperature on the atm grid
          strng = "Sea Ice Surface Temperature"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "K"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "Defined in the AGCM"
          ncstat = nf90_put_att( fid, vid, "comment", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("GTICE_ocn")
          !--- Sea Ice temperature on the ocn grid interpolated from GTICE_atm
          strng = "Sea Ice Surface Temperature"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "K"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "Interpolated from GTICE_atm"
          ncstat = nf90_put_att( fid, vid, "comment", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_TepIce_atm")
          !--- Ocean (LIM or CICE) Sea Ice Surface temperature on the atm grid
          !--- interpolated from O_TepIce
          strng = "Sea Ice Surface Temperature"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "K"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "Interpolated from O_TepIce"
          ncstat = nf90_put_att( fid, vid, "comment", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("GC_atm")
          !--- Ground cover mask
          strng = "Ground Cover Mask"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "-1 = land point, 0 = open ocean point, 1 = ice covered ocean point"
          ncstat = nf90_put_att( fid, vid, "comment", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("SIC_atm")
          !--- Sea ice thickness
          strng = "Sea Ice Thickness"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("SICN_atm")
          !--- Sea ice fraction
          strng = "Sea Ice Fraction"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("SNO_atm")
          !--- Snow thickness
          strng = "Snow Thickness"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("SNOW_atm")
          !--- Snow rate
          strng = "Total solid precipitation (snow rate)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m-2 s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("RES_atm")
          !--- Residual over lakes
          strng = "Residual Over Lakes"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("PMSL_atm")
          !--- Mean sea level pressure sent from the atm
          strng = "Mean Sea Level Pressure"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "hPa"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_MSLP")
          !--- Mean sea level pressure received by NEMO
          strng = "Mean Sea Level Pressure"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "atm"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("CO2_atm")
          !--- CO2 concentration sent from the atm
          strng = "CO2 Concentration (sent from AGCM)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "ppm"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_AtmCO2")
          !--- CO2 concentration received by NEMO
          strng = "CO2 Concentration (sent from AGCM)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "ppm"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("CO2flx_atm")
          !--- CO2 flux that is received by the AGCM
          strng = "CO2 Flux (sent from NEMO)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m2 s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_CO2FLX")
          !--- CO2 flux that is sent by NEMO
          strng = "CO2 Flux (sent from NEMO)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m2 s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("OUFS_atm")
          !--- Cell Average Zonal Atmosphere-Ocean Wind Stress
          strng = "Cell Average Zonal Atmospheric Wind Stress"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "N m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("UFSO_atm")
          !--- Zonal Atmosphere-Ocean Wind Stress
          strng = "Zonal Atmosphere-Ocean Wind Stress"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "N m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("UFSI_atm")
          !--- Zonal Atmosphere-Ice Wind Stress
          strng = "Zonal Atmosphere-Ice Wind Stress"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "N m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("OVFS_atm")
          !--- Cell Average Meridional Atmospheric Wind Stress
          strng = "Cell Average Meridional Atmospheric Wind Stress"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "N m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("VFSO_atm")
          !--- Meridional Atmosphere-Ocean Wind Stress
          strng = "Meridional Atmosphere-Ocean Wind Stress"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "N m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("VFSI_atm")
          !--- Meridional Atmosphere-Ice Wind Stress
          strng = "Meridional Atmosphere-Ice Wind Stress"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "N m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("OBEG_atm")
          !--- 
          strng = "Net Heat Flux"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("OBWG_atm")
          !--- Freshwater Budget (P-E)
          strng = "Freshwater Budget (P-E)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m-2 s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("ROFO_atm")
          !--- Runoff from land into ocean
          strng = "Land Surface Runoff"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m-2 s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("RIVO_atm")
          !--- Runoff from land into ocean
          strng = "Land Surface Runoff to Ocean"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m-2 s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("OFSG_atm")
          !--- Cell Average Net Solar Flux Over Ocean
          strng = "Cell Average Net Solar Flux Over Ocean"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("FSGO_atm")
          !--- Solar Flux Over Ocean
          strng = "Solar Flux Over Ocean"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("FSGI_atm")
          !--- Solar Flux Over Sea Ice
          strng = "Solar Flux Over Sea Ice"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("SWMX_atm")
          !--- 10 meter wind
          strng = "10 meter wind"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "m s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("BEGO_atm")
          !--- Total heat flux over ocean (positive down)
          strng = "Total heat flux over ocean"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("BEGI_atm")
          !--- Total heat flux over sea ice
          strng = "Total heat flux over sea ice"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("BWGO_atm")
          !--- Water Budget over ocean
          strng = "Ocean Freshwater Budget (P-E)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m-2 s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("BWGI_atm")
          !--- Water Budget over sea ice
          strng = "Sea Ice Freshwater Budget (P-E)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m-2 s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("RAIN_atm")
          !--- Rainfall Rate
          strng = "Total liquid precipitation"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m-2 s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("HFLI_atm")
          !--- Lsat * (Sublimation Over Sea Ice)
          strng = "Lsat * (Sublimation Over Sea Ice)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "Lsat = 2.834e6 J/kg  is the latent heat of sublimation for fresh ice or snow"
          ncstat = nf90_put_att( fid, vid, "comment", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("HSEA_atm")
          !--- Conductive heat flux over sea ice
          strng = "Conductive heat flux over sea ice"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("SLIM_atm")
          !--- Non-solar heat flux sensitivity to temperature over sea ice
          strng = "Non-solar sensitivity to temperature over sea ice (dQns/dT)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2 K-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("SLIMlw_atm")
          !--- Non-solar heat flux sensitivity to temperature over sea ice
          strng = "Longwave Term in dQns/dT"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2 K-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("SLIMsh_atm")
          !--- Non-solar heat flux sensitivity to temperature over sea ice
          strng = "Sensible Heat Term in dQns/dT"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2 K-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("SLIMlh_atm")
          !--- Non-solar heat flux sensitivity to temperature over sea ice
          strng = "Latent Heat Term in dQns/dT"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2 K-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_dQnsdT")
          !--- Non-solar heat flux sensitivity to temperature over sea ice
          strng = "Non-solar sensitivity to temperature over sea ice (dQns/dT)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2 K-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_dQlwdT")
          !--- Term in non-solar heat flux sensitivity to temperature over sea ice
          strng = "Longwave Term in dQns/dT"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2 K-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_dQshdT")
          !--- Term in non-solar heat flux sensitivity to temperature over sea ice
          strng = "Sensible Heat Term in dQns/dT"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2 K-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_dQlhdT")
          !--- Term in non-solar heat flux sensitivity to temperature over sea ice
          strng = "Latent Heat Term in dQns/dT"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2 K-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_OTaux1")
          !--- Atmosphere-Ocean wind stress
          strng = "Atmosphere-Ocean Wind Stress (x component)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "N m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_OTauy1")
          !--- Atmosphere-Ocean wind stress
          strng = "Atmosphere-Ocean Wind Stress (y component)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "N m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_ITaux1")
          !--- Atmosphere-Ice wind stress
          strng = "Atmosphere-Ice Wind Stress (x component)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "N m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_ITauy1")
          !--- Atmosphere-Ice wind stress
          strng = "Atmosphere-Ice Wind Stress (y component)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "N m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_QsrOce")
          !--- Solar heat flux over ocean
          strng = "Solar Heat Flux Over Ocean"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_QsrIce")
          !--- Solar heat flux over ice
          strng = "Solar Heat Flux Over Ice"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_QsrMix")
          !--- Mixed solar heat flux over ocean and ice
          strng = "Mixed Solar Heat Flux Over Ocean and Ice"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_QnsOce")
          !--- Non-Solar heat flux over ocean
          strng = "Non-Solar Heat Flux Over Ocean"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_QnsIce")
          !--- Non-Solar heat flux over ice
          strng = "Non-Solar Heat Flux Over Ice"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_QnsMix")
          !--- Mixed non-Solar heat flux over ocean and ice
          strng = "Mixed Non-Solar Heat Flux Over Ocean and Ice"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("OTotRain")
          !--- Total liquid precipitation
          strng = "Total Liquid Precipitation Over Ocean (rain)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m-2 s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("OTotSnow")
          !--- Total solid precipitation
          strng = "Solid Precipitation Over Ocean (snow)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m-2 s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("OTotEvap")
          !--- Total evaporation
          strng = "Total Evaporation Over Ocean"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m-2 s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("OIceEvap")
          !--- Evaporation over ice
          strng = "Solid Evaporation Over Ice (sublimation)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m-2 s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("OSubMPre")
          !--- Sublimation - liquid precip - solid precip
          strng = "Sublimation - Liquid precipitation - Solid precipitation"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m-2 s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("OISubMSn")
          !--- Ice freshwater budget
          strng = "Sea Ice freshwater budget (Sublimation - Snow)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m-2 s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("OOEvaMPr")
          !--- Ocean freshwater budget
          strng = "Ocean freshwater budget (Ocean Evaporation - Precipitation)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m-2 s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_Wind10")
          !--- 10 meter wind
          strng = "10 Meter Wind Speed Over Ocean"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "m s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_Runoff")
          !--- Runoff
          strng = "Runoff Into Ocean"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m-2 s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("OCalving")
          !--- Calving
          strng = "Calving of Ice into the Ocean"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "kg m-2 s-1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_TauMod")
          !--- Wind stress amplitude
          strng = "Wind Stress Amplitude Over Ocean"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "N m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("OTopMlt")
          !--- top melt
          strng = "Net Flux minus Conductive Flux at Top of Sea Ice"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("OBotMlt")
          !--- bottom melt
          strng = "Conductive Flux at Top of Sea Ice (positive down)"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "W m-2"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_SSTSST")
          !--- Ocean Surface Temperature
          strng = "Ocean Surface Temperature"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "K"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("O_TepIce")
          !--- Sea Ice Temperature
          strng = "Sea Ice Surface Temperature"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "K"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("OIceFrc")
          !--- Sea Ice Fraction
          strng = "Sea Ice Fraction"
          ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

          strng = "1"
          ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
          call nc_error_handler(fid, ncstat, 60, trim(strng2))

        case ("OIceTck")
          select case( trim(adjustl(sn_snd_thick_cldes)) )
            case ("weighted iwe and swe")
              !--- Sea Ice Water Equivalent
              strng = "Sea Ice Water Equivalent"
              ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
              call nc_error_handler(fid, ncstat, 60, trim(strng2))

              strng = "kg m-2"
              ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
              call nc_error_handler(fid, ncstat, 60, trim(strng2))

            case default
              !--- Sea Ice Thickness
              strng = "Sea Ice Thickness"
              ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
              call nc_error_handler(fid, ncstat, 60, trim(strng2))

              strng = "m"
              ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
              call nc_error_handler(fid, ncstat, 60, trim(strng2))
          end select

        case ("OSnwTck")
          select case( trim(adjustl(sn_snd_thick_cldes)) )
            case ("weighted iwe and swe")
              !--- Snow Water Equivalent over sea ice
              strng = "Snow Water Equivalent over Sea Ice"
              ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
              call nc_error_handler(fid, ncstat, 60, trim(strng2))

              strng = "kg m-2"
              ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
              call nc_error_handler(fid, ncstat, 60, trim(strng2))

            case default
              !--- Snow Thickness over Sea Ice
              strng = "Snow Thickness over Sea Ice"
              ncstat = nf90_put_att( fid, vid, "long_name", trim(strng) )
              call nc_error_handler(fid, ncstat, 60, trim(strng2))

              strng = "m"
              ncstat = nf90_put_att( fid, vid, "units", trim(strng) )
              call nc_error_handler(fid, ncstat, 60, trim(strng2))
          end select
          
      end select

    end subroutine nc_add_var_meta_data

!KEEP_FOR_NOW    This is more general than convert_corners but is not completely debugged
!KEEP_FOR_NOW    and is less efficient
!KEEP_FOR_NOW  
!KEEP_FOR_NOW    !***************************************************************************
!KEEP_FOR_NOW    !--- Convert 2D corner coordinate arrays dimensioned (4, grid_size)
!KEEP_FOR_NOW    !--- into a 1D Array without any repeating corner values 
!KEEP_FOR_NOW    !--- This assumes that there are exactly 4 corners per grid cell and all
!KEEP_FOR_NOW    !--- corners are ordered the same (counterclockwise) around each cell center
!KEEP_FOR_NOW    !***************************************************************************
!KEEP_FOR_NOW    subroutine convert_corner_arrays_to_1D(isSphere, dim1, dim2, &
!KEEP_FOR_NOW                                           cornerX2D, cornerY2D, cornerX, cornerY, rc)
!KEEP_FOR_NOW      implicit none
!KEEP_FOR_NOW
!KEEP_FOR_NOW      logical :: isSphere
!KEEP_FOR_NOW      integer(kind=4) :: dim1,dim2
!KEEP_FOR_NOW      real(kind=8) :: cornerX2D(:,:),cornerY2D(:,:)
!KEEP_FOR_NOW      real(kind=8) :: cornerX(:),cornerY(:)
!KEEP_FOR_NOW      integer :: rc
!KEEP_FOR_NOW
!KEEP_FOR_NOW      integer :: i,j
!KEEP_FOR_NOW      real(kind=8) :: tol=0.0000000001
!KEEP_FOR_NOW      integer :: topCorner
!KEEP_FOR_NOW      integer :: topRightCorner
!KEEP_FOR_NOW      integer :: BtmRightCorner
!KEEP_FOR_NOW      integer :: btmCorner
!KEEP_FOR_NOW      integer :: count,inPos,outPos
!KEEP_FOR_NOW      integer :: ip1,im1
!KEEP_FOR_NOW      logical :: matches
!KEEP_FOR_NOW     
!KEEP_FOR_NOW      ! make sure no dimensions are 0
!KEEP_FOR_NOW      if ((dim1 < 1) .or. (dim2 <1)) then
!KEEP_FOR_NOW        write(6,*)"convert_corner_arrays_to_1D: Currently can't handle a grid of width <1 in a dim."
!KEEP_FOR_NOW        call xit("convert_corner_arrays_to_1D",-1)
!KEEP_FOR_NOW      endif
!KEEP_FOR_NOW     
!KEEP_FOR_NOW      ! Handle 1 width cases
!KEEP_FOR_NOW      if ((dim1 == 1) .and. (dim2 == 1)) then
!KEEP_FOR_NOW        ! Put corner array into Grid in an order
!KEEP_FOR_NOW        ! so that it'll be recoverted to a 
!KEEP_FOR_NOW        ! Mesh element in the same order
!KEEP_FOR_NOW     
!KEEP_FOR_NOW         cornerX(1)=cornerX2D(1,1)
!KEEP_FOR_NOW         cornerY(1)=cornerY2D(1,1)    
!KEEP_FOR_NOW     
!KEEP_FOR_NOW         cornerX(2)=cornerX2D(2,1)
!KEEP_FOR_NOW         cornerY(2)=cornerY2D(2,1)    
!KEEP_FOR_NOW     
!KEEP_FOR_NOW         cornerX(3)=cornerX2D(4,1)
!KEEP_FOR_NOW         cornerY(3)=cornerY2D(4,1)    
!KEEP_FOR_NOW     
!KEEP_FOR_NOW         cornerX(4)=cornerX2D(3,1)
!KEEP_FOR_NOW         cornerY(4)=cornerY2D(3,1)    
!KEEP_FOR_NOW         return
!KEEP_FOR_NOW      endif
!KEEP_FOR_NOW     
!KEEP_FOR_NOW      if (dim1 == 1) then
!KEEP_FOR_NOW        write(6,*)"convert_corner_arrays_to_1D: Currently can't handle a grid thats width 1 in only 1st dim"
!KEEP_FOR_NOW        call xit("convert_corner_arrays_to_1D",-2)
!KEEP_FOR_NOW      endif
!KEEP_FOR_NOW     
!KEEP_FOR_NOW      if (dim2 == 1) then
!KEEP_FOR_NOW        write(6,*)"convert_corner_arrays_to_1D: Currently can't handle a grid thats width 1 in only 2nd dim"
!KEEP_FOR_NOW        call xit("convert_corner_arrays_to_1D",-3)
!KEEP_FOR_NOW      endif
!KEEP_FOR_NOW     
!KEEP_FOR_NOW     
!KEEP_FOR_NOW      ! Figure out which corner indice is the top row of corners
!KEEP_FOR_NOW      ! It won't match any of the neighbors corners
!KEEP_FOR_NOW      TopCorner=-1
!KEEP_FOR_NOW      do i=1,4
!KEEP_FOR_NOW     
!KEEP_FOR_NOW         ! See if it matches nbr to the right
!KEEP_FOR_NOW         matches=.false.
!KEEP_FOR_NOW         do j=1,4
!KEEP_FOR_NOW            if ((abs(cornerX2D(i,1)-cornerX2D(j,2))<tol) .and. &
!KEEP_FOR_NOW                (abs(cornerY2D(i,1)-cornerY2D(j,2))<tol)) then
!KEEP_FOR_NOW               matches=.true.
!KEEP_FOR_NOW               exit
!KEEP_FOR_NOW            endif
!KEEP_FOR_NOW         enddo
!KEEP_FOR_NOW         if (matches) cycle 
!KEEP_FOR_NOW     
!KEEP_FOR_NOW         ! See if it matches nbr to the below
!KEEP_FOR_NOW         matches=.false.
!KEEP_FOR_NOW         do j=1,4
!KEEP_FOR_NOW            if ((abs(cornerX2D(i,1)-cornerX2D(j,dim1+1))<tol) .and. &
!KEEP_FOR_NOW                (abs(cornerY2D(i,1)-cornerY2D(j,dim1+1))<tol)) then
!KEEP_FOR_NOW               matches=.true.
!KEEP_FOR_NOW               exit
!KEEP_FOR_NOW            endif
!KEEP_FOR_NOW         enddo
!KEEP_FOR_NOW         if (matches) cycle 
!KEEP_FOR_NOW     
!KEEP_FOR_NOW        ! See if it matches nbr to the below and to the right
!KEEP_FOR_NOW         matches=.false.
!KEEP_FOR_NOW         do j=1,4
!KEEP_FOR_NOW            if ((abs(cornerX2D(i,1)-cornerX2D(j,dim1+2))<tol) .and. &
!KEEP_FOR_NOW                (abs(cornerY2D(i,1)-cornerY2D(j,dim1+2))<tol)) then
!KEEP_FOR_NOW               matches=.true.
!KEEP_FOR_NOW               exit
!KEEP_FOR_NOW            endif
!KEEP_FOR_NOW         enddo
!KEEP_FOR_NOW         if (matches) cycle 
!KEEP_FOR_NOW     
!KEEP_FOR_NOW         ! Doesn't match anyone
!KEEP_FOR_NOW         TopCorner=i
!KEEP_FOR_NOW     
!KEEP_FOR_NOW         ! Exit the loop
!KEEP_FOR_NOW         exit
!KEEP_FOR_NOW      enddo
!KEEP_FOR_NOW     
!KEEP_FOR_NOW      ! Make sure we found a corner
!KEEP_FOR_NOW      if (TopCorner == -1) then
!KEEP_FOR_NOW        write(6,*)"convert_corner_arrays_to_1D: Bad corner array in SCRIP file"
!KEEP_FOR_NOW        call xit("convert_corner_arrays_to_1D",-4)
!KEEP_FOR_NOW      endif
!KEEP_FOR_NOW     
!KEEP_FOR_NOW     
!KEEP_FOR_NOW      ! Figure out which corner indice is the top right row of corners
!KEEP_FOR_NOW      ! It will match the top right, but not the bottom right
!KEEP_FOR_NOW      TopRightCorner=-1
!KEEP_FOR_NOW      do i=1,4
!KEEP_FOR_NOW     
!KEEP_FOR_NOW         ! See if it matches nbr to the right
!KEEP_FOR_NOW         matches=.false.
!KEEP_FOR_NOW         do j=1,4
!KEEP_FOR_NOW            if ((abs(cornerX2D(i,1)-cornerX2D(j,2))<tol) .and. &
!KEEP_FOR_NOW                (abs(cornerY2D(i,1)-cornerY2D(j,2))<tol)) then
!KEEP_FOR_NOW               matches=.true.
!KEEP_FOR_NOW               exit
!KEEP_FOR_NOW            endif
!KEEP_FOR_NOW         enddo
!KEEP_FOR_NOW         if (.not. matches) cycle 
!KEEP_FOR_NOW     
!KEEP_FOR_NOW         ! See if it matches nbr to the below right
!KEEP_FOR_NOW         matches=.false.
!KEEP_FOR_NOW         do j=1,4
!KEEP_FOR_NOW            if ((abs(cornerX2D(i,1)-cornerX2D(j,dim1+2))<tol) .and. &
!KEEP_FOR_NOW                (abs(cornerY2D(i,1)-cornerY2D(j,dim1+2))<tol)) then
!KEEP_FOR_NOW               matches=.true.
!KEEP_FOR_NOW               exit
!KEEP_FOR_NOW            endif
!KEEP_FOR_NOW         enddo
!KEEP_FOR_NOW     
!KEEP_FOR_NOW         ! correct matching so should be BtmCorner
!KEEP_FOR_NOW         if (.not. matches) then 
!KEEP_FOR_NOW            TopRightCorner=i
!KEEP_FOR_NOW            exit
!KEEP_FOR_NOW         endif
!KEEP_FOR_NOW      enddo
!KEEP_FOR_NOW     
!KEEP_FOR_NOW      ! Make sure we found a corner
!KEEP_FOR_NOW      if (TopRightCorner == -1) then
!KEEP_FOR_NOW        write(6,*)"convert_corner_arrays_to_1D: Bad corner array in SCRIP file"
!KEEP_FOR_NOW        call xit("convert_corner_arrays_to_1D",-5)
!KEEP_FOR_NOW      endif
!KEEP_FOR_NOW     
!KEEP_FOR_NOW      ! Figure out which corner indice is the bottom row of corners
!KEEP_FOR_NOW      ! It will match the one below , but not the one to the right
!KEEP_FOR_NOW      BtmCorner=-1
!KEEP_FOR_NOW      do i=1,4
!KEEP_FOR_NOW     
!KEEP_FOR_NOW         ! See if it matches nbr to the right
!KEEP_FOR_NOW         matches=.false.
!KEEP_FOR_NOW         do j=1,4
!KEEP_FOR_NOW            if ((abs(cornerX2D(i,1)-cornerX2D(j,2))<tol) .and. &
!KEEP_FOR_NOW                (abs(cornerY2D(i,1)-cornerY2D(j,2))<tol)) then
!KEEP_FOR_NOW               matches=.true.
!KEEP_FOR_NOW               exit
!KEEP_FOR_NOW            endif
!KEEP_FOR_NOW         enddo
!KEEP_FOR_NOW         if (matches) cycle 
!KEEP_FOR_NOW     
!KEEP_FOR_NOW         ! See if it matches nbr to the below
!KEEP_FOR_NOW         matches=.false.
!KEEP_FOR_NOW         do j=1,4
!KEEP_FOR_NOW            if ((abs(cornerX2D(i,1)-cornerX2D(j,dim1+1))<tol) .and. &
!KEEP_FOR_NOW                (abs(cornerY2D(i,1)-cornerY2D(j,dim1+1))<tol)) then
!KEEP_FOR_NOW               matches=.true.
!KEEP_FOR_NOW               exit
!KEEP_FOR_NOW            endif
!KEEP_FOR_NOW         enddo
!KEEP_FOR_NOW     
!KEEP_FOR_NOW         ! correct matching so should be BtmCorner
!KEEP_FOR_NOW         if (matches) then 
!KEEP_FOR_NOW            BtmCorner=i
!KEEP_FOR_NOW            exit
!KEEP_FOR_NOW         endif
!KEEP_FOR_NOW     
!KEEP_FOR_NOW      enddo
!KEEP_FOR_NOW     
!KEEP_FOR_NOW      ! Make sure we found a corner
!KEEP_FOR_NOW      if (BtmCorner == -1) then
!KEEP_FOR_NOW        write(6,*)"convert_corner_arrays_to_1D: Bad corner array in SCRIP file"
!KEEP_FOR_NOW        call xit("convert_corner_arrays_to_1D",-6)
!KEEP_FOR_NOW      endif
!KEEP_FOR_NOW     
!KEEP_FOR_NOW     
!KEEP_FOR_NOW      ! Figure out which corner indice is the bottom right row of corners
!KEEP_FOR_NOW      ! It will match the bottom right, but not the top right
!KEEP_FOR_NOW      BtmRightCorner=-1
!KEEP_FOR_NOW      do i=1,4
!KEEP_FOR_NOW     
!KEEP_FOR_NOW         ! eliminate all other possibilities
!KEEP_FOR_NOW         if (i == TopCorner) cycle
!KEEP_FOR_NOW         if (i == TopRightCorner) cycle
!KEEP_FOR_NOW         if (i == BtmCorner) cycle
!KEEP_FOR_NOW     
!KEEP_FOR_NOW         BtmRightCorner=i
!KEEP_FOR_NOW      enddo
!KEEP_FOR_NOW     
!KEEP_FOR_NOW      ! Make sure we found a corner
!KEEP_FOR_NOW      if (BtmRightCorner == -1) then
!KEEP_FOR_NOW        write(6,*)"convert_corner_arrays_to_1D: Bad corner array in SCRIP file"
!KEEP_FOR_NOW        call xit("convert_corner_arrays_to_1D",-7)
!KEEP_FOR_NOW      endif
!KEEP_FOR_NOW     
!KEEP_FOR_NOW     
!KEEP_FOR_NOWif ( .false. ) then ! DBG1
!KEEP_FOR_NOW
!KEEP_FOR_NOW      ! Error check corner info from file to make sure corners are consistent throughout file
!KEEP_FOR_NOW      if (isSphere) then
!KEEP_FOR_NOW        write(*,*) 
!KEEP_FOR_NOW        write(*,*) "Error checking spherical Grid..."
!KEEP_FOR_NOW        
!KEEP_FOR_NOW        ! Init to starting pos
!KEEP_FOR_NOW        inPos=0
!KEEP_FOR_NOW        
!KEEP_FOR_NOW        ! Check rows up to btm row
!KEEP_FOR_NOW        do i=1,dim2-1
!KEEP_FOR_NOW           do j=1,dim1-1
!KEEP_FOR_NOW              inPos=inPos+1
!KEEP_FOR_NOW              if ((cornerX2D(TopRightCorner,inPos) .ne. cornerX2D(TopCorner,inPos+1)) .and. &
!KEEP_FOR_NOW                   (cornerY2D(TopRightCorner,inPos) .ne. cornerY2D(TopCorner,inPos+1))) then
!KEEP_FOR_NOW                 write(*,*) "TopRightCorner of ",inPos," doesn't match TopCorner of ",inPos+1
!KEEP_FOR_NOW                 write(*,*) cornerX2D(TopRightCorner,inPos),cornerY2D(TopRightCorner,inPos)," .ne. ", &
!KEEP_FOR_NOW                            cornerX2D(TopCorner,inPos+1),cornerY2D(TopCorner,inPos+1)
!KEEP_FOR_NOW                 !              stop
!KEEP_FOR_NOW              endif
!KEEP_FOR_NOW              
!KEEP_FOR_NOW              if ((cornerX2D(BtmRightCorner,inPos) .ne. cornerX2D(BtmCorner,inPos+1)) .and. &
!KEEP_FOR_NOW                   (cornerY2D(BtmRightCorner,inPos) .ne. cornerY2D(BtmCorner,inPos+1))) then
!KEEP_FOR_NOW                 write(*,*) "BtmRightCorner of ",inPos," doesn't match BtmCorner of ",inPos+1
!KEEP_FOR_NOW                 !           stop
!KEEP_FOR_NOW              endif
!KEEP_FOR_NOW              
!KEEP_FOR_NOW              if ((cornerX2D(BtmCorner,inPos) .ne. cornerX2D(TopCorner,inPos+dim1)) .and. &
!KEEP_FOR_NOW                   (cornerY2D(BtmCorner,inPos) .ne. cornerY2D(TopCorner,inPos+dim1))) then
!KEEP_FOR_NOW                 write(*,*) "BtmCorner of ",inPos," doesn't match TopCorner of ",inPos+dim1
!KEEP_FOR_NOW                 !           stop
!KEEP_FOR_NOW              endif
!KEEP_FOR_NOW              
!KEEP_FOR_NOW              if ((cornerX2D(BtmRightCorner,inPos) .ne. cornerX2D(TopRightCorner,inPos+dim1)) .and. &
!KEEP_FOR_NOW                   (cornerY2D(BtmRightCorner,inPos) .ne. cornerY2D(TopRightCorner,inPos+dim1))) then
!KEEP_FOR_NOW                 write(*,*) "BtmRightCorner of ",inPos," doesn't match TopRightCorner of ",inPos+dim1
!KEEP_FOR_NOW                 !           stop
!KEEP_FOR_NOW              endif
!KEEP_FOR_NOW           enddo
!KEEP_FOR_NOW           
!KEEP_FOR_NOW           ! Check last point in row with beginning of row 
!KEEP_FOR_NOW           inPos=inPos+1
!KEEP_FOR_NOW           if ((cornerX2D(TopRightCorner,inPos) .ne. cornerX2D(TopCorner,inPos-dim1+1)) .and. &
!KEEP_FOR_NOW                (cornerY2D(TopRightCorner,inPos) .ne. cornerY2D(TopCorner,inPos-dim1+1))) then
!KEEP_FOR_NOW              write(*,*) "TopRightCorner of ",inPos," doesn't match TopCorner of ",inPos-dim1+1
!KEEP_FOR_NOW              !        stop
!KEEP_FOR_NOW           endif
!KEEP_FOR_NOW           
!KEEP_FOR_NOW           if ((cornerX2D(BtmRightCorner,inPos) .ne. cornerX2D(BtmCorner,inPos-dim1+1)) .and. &
!KEEP_FOR_NOW                (cornerY2D(BtmRightCorner,inPos) .ne. cornerY2D(BtmCorner,inPos-dim1+1))) then
!KEEP_FOR_NOW              write(*,*) "BtmRightCorner of ",inPos," doesn't match BtmCorner of ",inPos-dim1+1
!KEEP_FOR_NOW              stop
!KEEP_FOR_NOW           endif
!KEEP_FOR_NOW           
!KEEP_FOR_NOW           
!KEEP_FOR_NOW           if ((cornerX2D(BtmCorner,inPos) .ne. cornerX2D(TopCorner,inPos+dim1)) .and. &
!KEEP_FOR_NOW                (cornerY2D(BtmCorner,inPos) .ne. cornerY2D(TopCorner,inPos+dim1))) then
!KEEP_FOR_NOW              write(*,*) "BtmCorner of ",inPos," doesn't match TopCorner of ",inPos+dim1
!KEEP_FOR_NOW              !        stop
!KEEP_FOR_NOW           endif
!KEEP_FOR_NOW           
!KEEP_FOR_NOW           if ((cornerX2D(BtmRightCorner,inPos) .ne. cornerX2D(TopRightCorner,inPos+dim1)) .and. &
!KEEP_FOR_NOW                (cornerY2D(BtmRightCorner,inPos) .ne. cornerY2D(TopRightCorner,inPos+dim1))) then
!KEEP_FOR_NOW              write(*,*) "BtmRightCorner of ",inPos," doesn't match TopRightCorner of ",inPos+dim1
!KEEP_FOR_NOW              !        stop
!KEEP_FOR_NOW           endif
!KEEP_FOR_NOW        enddo
!KEEP_FOR_NOW     
!KEEP_FOR_NOW     
!KEEP_FOR_NOW        ! Check bottom row
!KEEP_FOR_NOW        do j=1,dim1-1
!KEEP_FOR_NOW           inPos=inPos+1
!KEEP_FOR_NOW           if ((cornerX2D(TopRightCorner,inPos) .ne. cornerX2D(TopCorner,inPos+1)) .and. &
!KEEP_FOR_NOW                (cornerY2D(TopRightCorner,inPos) .ne. cornerY2D(TopCorner,inPos+1))) then
!KEEP_FOR_NOW              write(*,*) "TopRightCorner of ",inPos," doesn't match TopCorner of ",inPos+1
!KEEP_FOR_NOW              !              stop
!KEEP_FOR_NOW           endif
!KEEP_FOR_NOW           
!KEEP_FOR_NOW           if ((cornerX2D(BtmRightCorner,inPos) .ne. cornerX2D(BtmCorner,inPos+1)) .and. &
!KEEP_FOR_NOW                (cornerY2D(BtmRightCorner,inPos) .ne. cornerY2D(BtmCorner,inPos+1))) then
!KEEP_FOR_NOW              write(*,*) "BtmRightCorner of ",inPos," doesn't match BtmCorner of ",inPos+1
!KEEP_FOR_NOW              !           stop
!KEEP_FOR_NOW           endif
!KEEP_FOR_NOW        enddo
!KEEP_FOR_NOW     
!KEEP_FOR_NOW     
!KEEP_FOR_NOW        ! Check last point in row with beginning of row 
!KEEP_FOR_NOW        inPos=inPos+1
!KEEP_FOR_NOW        if ((cornerX2D(TopRightCorner,inPos) .ne. cornerX2D(TopCorner,inPos-dim1+1)) .and. &
!KEEP_FOR_NOW             (cornerY2D(TopRightCorner,inPos) .ne. cornerY2D(TopCorner,inPos-dim1+1))) then
!KEEP_FOR_NOW           write(*,*) "TopRightCorner of ",inPos," doesn't match TopCorner of ",inPos-dim1+1
!KEEP_FOR_NOW           !        stop
!KEEP_FOR_NOW        endif
!KEEP_FOR_NOW        
!KEEP_FOR_NOW        if ((cornerX2D(BtmRightCorner,inPos) .ne. cornerX2D(BtmCorner,inPos-dim1+1)) .and. &
!KEEP_FOR_NOW             (cornerY2D(BtmRightCorner,inPos) .ne. cornerY2D(BtmCorner,inPos-dim1+1))) then
!KEEP_FOR_NOW           write(*,*) "BtmRightCorner of ",inPos," doesn't match BtmCorner of ",inPos-dim1+1
!KEEP_FOR_NOW           !     stop
!KEEP_FOR_NOW        endif
!KEEP_FOR_NOW     else
!KEEP_FOR_NOW        ! TODO: Check regional grid
!KEEP_FOR_NOW           
!KEEP_FOR_NOW           
!KEEP_FOR_NOW     endif
!KEEP_FOR_NOW
!KEEP_FOR_NOWendif ! DBG1
!KEEP_FOR_NOW     
!KEEP_FOR_NOW      ! Set Corner info
!KEEP_FOR_NOW      if (isSphere) then
!KEEP_FOR_NOW         count=size(cornerX2D,2)
!KEEP_FOR_NOW         do i=1,count
!KEEP_FOR_NOW            cornerX(i)=cornerX2D(TopCorner,i)
!KEEP_FOR_NOW            cornerY(i)=cornerY2D(TopCorner,i)
!KEEP_FOR_NOW         enddo
!KEEP_FOR_NOW          
!KEEP_FOR_NOW         do i=1,dim1
!KEEP_FOR_NOW            cornerX(i+count)=cornerX2D(BtmCorner,count-dim1+i)
!KEEP_FOR_NOW            cornerY(i+count)=cornerY2D(BtmCorner,count-dim1+i)
!KEEP_FOR_NOW         enddo
!KEEP_FOR_NOW      else
!KEEP_FOR_NOW         ! Set Corner info
!KEEP_FOR_NOW         inPos=0
!KEEP_FOR_NOW         outPos=0
!KEEP_FOR_NOW         do i=1,dim2
!KEEP_FOR_NOW            do j=1,dim1
!KEEP_FOR_NOW               inPos=inPos+1
!KEEP_FOR_NOW               outPos=outPos+1
!KEEP_FOR_NOW               cornerX(outPos)=cornerX2D(TopCorner,inPos)
!KEEP_FOR_NOW               cornerY(outPos)=cornerY2D(TopCorner,inPos)
!KEEP_FOR_NOW            enddo
!KEEP_FOR_NOW            outPos=outPos+1
!KEEP_FOR_NOW            cornerX(outPos)=cornerX2D(TopRightCorner,inPos)
!KEEP_FOR_NOW            cornerY(outPos)=cornerY2D(TopRightCorner,inPos)
!KEEP_FOR_NOW         enddo
!KEEP_FOR_NOW          
!KEEP_FOR_NOW         inPos=inPos-dim1
!KEEP_FOR_NOW         do i=1,dim1
!KEEP_FOR_NOW            inPos=inPos+1
!KEEP_FOR_NOW            outPos=outPos+1
!KEEP_FOR_NOW            cornerX(outPos)=cornerX2D(BtmCorner,inPos)
!KEEP_FOR_NOW            cornerY(outPos)=cornerY2D(BtmCorner,inPos)
!KEEP_FOR_NOW         enddo
!KEEP_FOR_NOW          
!KEEP_FOR_NOW         outPos=outPos+1
!KEEP_FOR_NOW         cornerX(outPos)=cornerX2D(BtmRightCorner,inPos)
!KEEP_FOR_NOW         cornerY(outPos)=cornerY2D(BtmRightCorner,inPos)
!KEEP_FOR_NOW      endif
!KEEP_FOR_NOW       
!KEEP_FOR_NOW      ! return success
!KEEP_FOR_NOW      rc=1
!KEEP_FOR_NOW     
!KEEP_FOR_NOW    end subroutine convert_corner_arrays_to_1D

end module ccc_netcdf
