module regrid_mod

  use ccc_esmf
  use grid_desc

  public

  !--- used to store command line arguments
  integer*4 iargc,iarg
  integer nargs
  character(256) :: clarg, pname, pval
  character(256), allocatable :: clargs(:)

  integer :: idx1, idx2, verbose=1
  integer :: nrec, kfile, nfiles=0
  integer :: nx, ny, nwrds, nx1, ny1, nx2, ny2

  character(256) :: curr_file, files(100)
  character(256) :: fout, forig
  logical :: usr_def_fout
  character(256) :: mask_gg, ocn_desc, nemo_coords

  character(len=32) :: overwrite_with = " "

  logical :: use_timers = .true.

  integer :: start_ibuf2, stop_ibuf2

  !--- Flag the interpolastion of a test function rather than user supplied data
  logical :: check_remap = .false.

  !--- Flags the creation of grid description files only
  logical :: grid_desc_only = .false.

  !--- Mask or don't mask interpolation
  logical :: src_masked=.true.
  logical :: dst_masked=.true.
  integer :: curr_mask_type = 0

  logical :: write_orig_fld = .false.
  logical :: write_reinterp_diff = .false.
  logical :: add_overlap_lon = .false.

  !--- olaplon determine if the grid description files created here will
  !--- contain overlap longitudes (atm and old ocn grids only)
  !--- NOTE: grids retrieved from these files may or may not contain the overlap lon
  logical :: olaplon = .true.

  !--- The maximum number of records that may be read from a single input file
  integer :: max_rec_read = 9999999

  integer :: nvlist=0
  character(256) :: vlist(100)

  !--- Netcdf related
  character(128) :: nc_file, nc_name
  integer(kind=8) :: curr_time
  logical :: nc_out = .false.

  character(64) :: method = "regrid"
  integer(ESMF_KIND_I) :: lnum
  integer(ESMF_KIND_I) :: rc

  !--- Name/alias of src and dst grid and grid description files
  character(256) :: src_grid, dst_grid

  !--- The type of interpolation to be done
  character(64) :: regrid_method = 'bilinear'

  !--- The type of normalization to be done for a conservative remap
  character(64) :: this_norm_opt = 'fracarea'

  real(kind=8), pointer :: src_data(:)
  real(kind=8), pointer :: dst_data(:)

  integer(kind=8), pointer :: itmp_data(:)
  real(kind=8), pointer    :: rtmp_data(:)

  !--- For grid definition calls
  integer :: id_src_grid, id_dst_grid
  logical :: use_src_mask, use_dst_mask

  integer :: src_grid_dims(2), dst_grid_dims(2)
  integer :: src_grid_edges(2,2), dst_grid_edges(2,2)
  integer :: src_grid_size, dst_grid_size
  integer(ESMF_KIND_I4), pointer :: src_grid_mask(:,:)
  integer(ESMF_KIND_I4), pointer :: dst_grid_mask(:,:)
  real(ESMF_KIND_R8), pointer :: src_grid_area(:,:)
  real(ESMF_KIND_R8), pointer :: dst_grid_area(:,:)

  real(ESMF_KIND_R8), pointer :: src_cent_lon(:,:)
  real(ESMF_KIND_R8), pointer :: src_cent_lat(:,:)
  real(ESMF_KIND_R8), pointer :: dst_cent_lon(:,:)
  real(ESMF_KIND_R8), pointer :: dst_cent_lat(:,:)
  real(ESMF_KIND_R8), pointer :: src_corn_lon(:,:)
  real(ESMF_KIND_R8), pointer :: src_corn_lat(:,:)
  real(ESMF_KIND_R8), pointer :: dst_corn_lon(:,:)
  real(ESMF_KIND_R8), pointer :: dst_corn_lat(:,:)

  !--- Some constants of default ESMF integer type
  integer(ESMF_KIND_I), parameter :: zero=0
  integer(ESMF_KIND_I), parameter ::  one=1
  integer(ESMF_KIND_I), parameter ::  two=2

  character(512) :: strng

  integer :: regrid_id

  real(ESMF_KIND_R8), pointer, dimension(:,:) :: src_area_data => null()
  real(ESMF_KIND_R8), pointer, dimension(:,:) :: dst_area_data => null()

  logical :: create_frac_land_mask = .false.
  logical :: frac_land_mask_add_olap = .false.
  real(kind=8), pointer :: src_fld(:), dst_fld(:)

  real(kind=8), pointer :: src_frac(:), dst_frac(:), src_area(:), dst_area(:)
  integer :: nx_sfrac, ny_sfrac, nx_dfrac, ny_dfrac
  integer :: nx_sarea, ny_sarea, nx_darea, ny_darea

  type(cpl_vinfo_t), pointer :: src_vinfo, dst_vinfo

  integer, pointer :: src_mask_values(:), dst_mask_values(:)

  !--- Hard code the min/max values used in test_function
  real(kind=8), parameter :: test_func_xmin=0.0_8, test_func_xmax=6.0_8
  real(kind=8), parameter :: test_func_ymin=0.0_8, test_func_ymax=6.0_8

  !--- Used by test_function
  logical, save :: interp_mask=.false., interp_sicn=.false.

  !--- Buffer and work space used by CCCma io routines
  integer, parameter :: maxx=8392704  !---2049x4096
  real               :: wrk(maxx)
  integer            :: ibuf(8)
  integer            :: idat(2*maxx)
  common /icom/ ibuf, idat

contains

  subroutine process_clg()
    !--- Determine the number of command line arguments
    nargs=iargc()

    !--- If there were no command line args then print
    !--- complete usage info to stdout and quit
    if (nargs.eq.0) call usage(0)

    if (nargs.lt.0) then
      !--- This should never happen, but ...
      write(6,'(a,i8)')'regrid: nargs is out of range. nargs = ',nargs
      call usage(-1)
    endif

    if (nargs.gt.0) then
      allocate ( clargs(nargs) )

      !--- Populate clargs with all command line arguments
      do iarg=1,nargs
        call getarg(iarg,clargs(iarg))
      enddo

      !--- Process command line arguments
      do idx1=1,nargs
        clarg=' '
        clarg=clargs(idx1)
        if (clarg(1:1).eq.'-') then
          !--- This is an option
          if (len_trim(clarg).eq.2 .and. clarg(1:2).eq.'-v') then
            verbose = verbose + 1
          else if (len_trim(clarg).eq.2 .and. clarg(1:2).eq.'-c') then
            grid_desc_only = .true.
          else if (len_trim(clarg).eq.3 .and. clarg(1:3).eq.'-wo') then
            write_orig_fld = .true.
          else if (len_trim(clarg).eq.3 .and. clarg(1:3).eq.'-w2') then
            write_reinterp_diff = .true.
          else if (len_trim(clarg).eq.3 .and. clarg(1:3).eq.'-go') then
            add_overlap_lon = .true.
          else if (len_trim(clarg).eq.6 .and. clarg(1:6).eq.'-nolap') then
            olaplon = .false.
          else if (len_trim(clarg).eq.3 .and. clarg(1:3).eq.'-nc') then
            nc_out = .true.
          else if (len_trim(clarg).eq.11 .and. clarg(1:11).eq.'-nosrc_mask') then
            !--- Do not use any mask on the source grid during interpolation
            src_masked = .false.
          else if (len_trim(clarg).eq.11 .and. clarg(1:11).eq.'-nodst_mask') then
            !--- Do not use any mask on the destination grid during interpolation
            dst_masked = .false.
          else if (len_trim(clarg).eq.6 .and. clarg(1:6).eq.'-check') then
            !--- Interpolate a test function rather than the data supplied
            check_remap = .true.
          else if (len_trim(clarg).eq.12 .and. clarg(1:12).eq.'-interp_mask') then
            !--- Interpolate an example land mask in the subroutine test_grid
            interp_mask = .true.
          else if (len_trim(clarg).eq.12 .and. clarg(1:12).eq.'-interp_sicn') then
            !--- Interpolate an example SICN in the subroutine test_grid
            interp_sicn = .true.
          else if (len_trim(clarg).eq.7 .and. clarg(1:7).eq.'-flmask') then
            !--- Create the AGCM fractional land mask from the ORCA1 land mask
            !--- and add an ovelap longitude to the output file
            create_frac_land_mask   = .true.
            frac_land_mask_add_olap = .true.
          else if (len_trim(clarg).eq.15 .and. clarg(1:15).eq.'-flmask_no_olap') then
            !--- Create the AGCM fractional land mask from the ORCA1 land mask
            !--- but do not add an ovelap longitude to the output file
            create_frac_land_mask   = .true.
            frac_land_mask_add_olap = .false.
          else
            write(6,*)'regrid: Unknown command line option ',trim(clarg)
            call usage(-2)
          endif
  
        else if (index(clarg,'=').gt.0) then
          !--- This argument is of the form vname=value
          idx2=index(clarg,'=')
          pname = " "
          pval  = " "
          pname = trim(clarg(1:idx2-1))
          pval  = trim(clarg(idx2+1:))
          if (len_trim(pval).le.0) then
            write(6,'(3a)')'regrid: Command line parameter ',trim(pname),' requres a value.'
            call usage(-3)
          endif

          select case (trim(adjustl(pname)))
            case ('verbose')
              !--- Set verbosity level
              pval=trim(adjustl(pval))
              if (verify(trim(pval),'0123456789') .ne. 0) then
                !--- Verify that pval is an integer
                write(6,*)'regrid: ',trim(pval),' is not an integer.'
                call usage(-3)
              endif
              read(pval,*,iostat=ios)verbose
              if (ios.ne.0) then
                write(6,*)'regrid: Error ',ios,' reading internal file ',trim(pval)
                call usage(-3)
              endif

            case ('max_read')
              !--- Set the maximum number of records that may be read from a single input file
              pval=trim(adjustl(pval))
              if (verify(trim(pval),'0123456789') .ne. 0) then
                !--- Verify that pval is an integer
                write(6,*)'regrid: ',trim(pval),' is not an integer.'
                call usage(-3)
              endif
              read(pval,*,iostat=ios)max_rec_read
              if (ios.ne.0) then
                write(6,*)'regrid: Error ',ios,' reading internal file ',trim(pval)
                call usage(-3)
              endif

            case ('start')
              !--- Set the starting value for ibuf2 in records read from the input file
              pval=trim(adjustl(pval))
              if (verify(trim(pval),'0123456789') .ne. 0) then
                !--- Verify that pval is an integer
                write(6,*)'regrid: ',trim(pval),' is not an integer.'
                call usage(-3)
              endif
              read(pval,*,iostat=ios)start_ibuf2
              if (ios.ne.0) then
                write(6,*)'regrid: Error ',ios,' reading internal file ',trim(pval)
                call usage(-3)
              endif

            case ('stop')
              !--- Set the stopping value for ibuf2 in records read from the input file
              pval=trim(adjustl(pval))
              if (verify(trim(pval),'0123456789') .ne. 0) then
                !--- Verify that pval is an integer
                write(6,*)'regrid: ',trim(pval),' is not an integer.'
                call usage(-3)
              endif
              read(pval,*,iostat=ios)stop_ibuf2
              if (ios.ne.0) then
                write(6,*)'regrid: Error ',ios,' reading internal file ',trim(pval)
                call usage(-3)
              endif

            case ('var', 'vars')
              !--- This should be a white space separated list of CCCma variable names
              !--- Any name longer than 4 chars will be silently truncated
              call split_on_delim(vlist,nvlist,trim(adjustl(pval))," ")

            case ('src_grid', 'src')
              !--- The name of file containing a description of the source grid
              !--- or a valid grid id (see below for a description of valid grid ids)
              src_grid = " "
              src_grid = trim(adjustl(pval))

            case ('dst_grid', 'dst')
              !--- The name of file containing a description of the destination grid
              !--- or a valid grid id (see below for a description of valid grid ids)
              dst_grid = " "
              dst_grid = trim(adjustl(pval))

            case ('ocn_desc')
              !--- The name of a CCCma binary file containing ocean data description info
              !--- This file must define: LON, LONH, LAT, LATH, Z, ZH, ZBOT
              ocn_desc = " "
              ocn_desc = trim(adjustl(pval))

            case ('mask_gg')
              !--- The name of a CCCma binary file containing a gaussian grid land mask
              !---   -1 for land, anything else is not-land
              !--- This file must be consistent with the gaussian grid to be created
              mask_gg = " "
              mask_gg = trim(adjustl(pval))

            case ('nemo_coords')
              !--- nemo_coords is the name of an existing NEMO netcdf file containing
              !--- ORCA coordinate values and land mask in variables named as follows
              !---     glamt, gphit, tmask :: T grid coords and land mask
              !---     glamu, gphiu, umask :: U grid coords and land mask
              !---     glamv, gphiv, vmask :: V grid coords and land mask
              !---     glamf, gphif, fmask :: F grid coords and land mask
              !--- The glam* variables are longitudes
              !--- The gphi* variables are latitudes
              !--- The *mask variable is optional and will default to mask == 1
              !--- if not present in the file
              nemo_coords = " "
              nemo_coords = trim(adjustl(pval))

            case ('map_type', 'method' )
              !--- The type of interpolation to be done
              regrid_method = " "
              regrid_method = lowerc( trim(adjustl(pval)) )

            case ('norm_opt')
              !--- The type of normalization to be done for a conservative mapping
              !--- Valid values are: fracarea, destarea, none (case is ignored)
              this_norm_opt = " "
              this_norm_opt = lowerc( trim(adjustl(pval)) )

            case ('overwrite')
              !--- Identify a field that will overwrite the output array from new_map
              !--- Valid values are: grid_frac, grid_area
              overwrite_with = " "
              overwrite_with = lowerc( trim(adjustl(pval)) )
              select case ( trim(adjustl(overwrite_with)) )
                case ("grid_frac")
                case ("grid_area")
                case default
                  write(6,*)'regrid: Invalid value for overwrite = ',trim(overwrite_with)
                  call usage(-3)
              end select

            case ('out')
              !--- The output file name
              fout = " "
              fout = trim(adjustl(pval))
              usr_def_fout = .true.

            case default
              !--- Unknown variable definition
              write(6,*)'regrid: Invalid command line arg ',trim(clarg)
              call usage(-3)

          end select

        else
          !--- This is assumed to be a file name
          nfiles = nfiles+1
          files(nfiles) = trim(adjustl(clarg))

        endif
      enddo

      deallocate ( clargs )
    endif

  end subroutine process_clg

  subroutine usage(mode)

    implicit none

    integer, intent(in) :: mode

    !--- local
    character(120) :: lines(500), strng
    integer        :: k,nlines,nlines_out

    !--- Define a usage message
    nlines=1
    lines(nlines) = 'Usage: regrid [options|defs] src_grid=FILE dst_grid=FILE data_file ...'
    nlines=nlines+1
    lines(nlines) = '  Interpolate input data onto another grid.'
    nlines=nlines+1
    lines(nlines) = ' '
    nlines=nlines+1
    lines(nlines) = 'Options:'
    nlines=nlines+1
    lines(nlines) = '  -nosrc_mask ..Do not mask source grid points'
    nlines=nlines+1
    lines(nlines) = '  -nodst_mask ..Do not mask destination grid points'
    nlines=nlines+1
    lines(nlines) = '               By default the mask supplied in the input grid description files is used'
    nlines=nlines+1
    lines(nlines) = '  -nolap     ..grid description files created here will not contain overlap longitude'
    nlines=nlines+1
    lines(nlines) = '  -wo        ..Write original fields to a second output file'
    nlines=nlines+1
    lines(nlines) = '  -w2        ..Write a second field to the output file containing the difference'
    nlines=nlines+1
    lines(nlines) = '               between the original field and a twice interpolated field (onto'
    nlines=nlines+1
    lines(nlines) = '               the target grid then back onto the original grid)'
    nlines=nlines+1
    lines(nlines) = '  -c         ..Create the src grid description file and exit'
    nlines=nlines+1
    lines(nlines) = '  -flmask     ..Create a CCCma binary file containing a fractional land mask, then exit'
    nlines=nlines+1
    lines(nlines) = '  -flmask_no_olap ..Same as above but without any cyclic longitude'
    nlines=nlines+1
    lines(nlines) = '               By default the 2 grids involved are T63 gaussian and ORCA1 but this can'
    nlines=nlines+1
    lines(nlines) = '               be changed by adding either src=FILE or dst=FILE to the command line.'
    nlines=nlines+1
    lines(nlines) = '  -check     ..Interpolate a test function rather than the supplied data'
    nlines=nlines+1
    lines(nlines) = '               The output will be the error in the interpolation'
    nlines=nlines+1
    lines(nlines) = '  -nc        ..Output is to a netcdf file (default is to a CCCma binary file)'
    nlines=nlines+1
    lines(nlines) = '  -v         ..Increment verbosity'
    nlines=nlines+1
    lines(nlines) = ' '
    nlines=nlines+1
    lines(nlines) = 'Definitions:'
    nlines=nlines+1
    lines(nlines) = '  src=FILE         ..A netcdf file containing a description of the source grid'
    nlines=nlines+1
    lines(nlines) = '  dst=FILE         ..A netcdf file containing a description of the destination grid'
    nlines=nlines+1
    lines(nlines) = '                     src and dst may also contain a grid ID, in which'
    nlines=nlines+1
    lines(nlines) = '                     case the grid description file will be created at run time.'
    nlines=nlines+1
    lines(nlines) = '                     Valid IDs: t31, t47, t63, t95, t127, orca1, ocnt_256x192, ocnv_256x192'
    nlines=nlines+1
    lines(nlines) = '  mask_gg=FILE     ..A CCCma binary file containing a gaussian grid land mask'
    nlines=nlines+1
    lines(nlines) = '                     The mask will contain -1 for land, anything else is not-land'
    nlines=nlines+1
    lines(nlines) = '  ocn_desc=FILE    ..A CCCma binary file containing ocean data description info'
    nlines=nlines+1
    lines(nlines) = '                     This file must define: LON, LONH, LAT, LATH, Z, ZH, ZBOT'
    nlines=nlines+1
    lines(nlines) = '  nemo_coords=FILE ..A netcdf file containing a NEMO ORCA grid description'
    nlines=nlines+1
    lines(nlines) = '                     This file must define: glam_, gphi_, _mask where _ = t,u,v,f'
    nlines=nlines+1
    lines(nlines) = '  method=STRING    ..The type of interpolation to be done. Valid values are:'
    nlines=nlines+1
    lines(nlines) = '                     conservative|cv, bilinear|bl, patch, nearest_stod|stod, nearest_dtos|dtos'
    nlines=nlines+1
    lines(nlines) = '  norm_opt=STRING  ..The type of normalization to be done for a conservative remap'
    nlines=nlines+1
    lines(nlines) = '                     Valid values are: fracarea, destarea, none'
    nlines=nlines+1
    lines(nlines) = '  var[s]=LIST      ..A whitespace separated list of CCCma variable names (ibuf3 values)'
    nlines=nlines+1
    lines(nlines) = '                     to extract from input files. All other records are ignored.'
    nlines=nlines+1
    lines(nlines) = '  start=INT        ..Records with ibuf2 values less than start will be ignored'
    nlines=nlines+1
    lines(nlines) = '  stop=INT         ..Records with ibuf2 values greater than stop will be ignored'

    nlines_out = 0

    if (mode.ge.0) nlines_out = nlines

    do k=1,nlines_out
      write(6,*)trim(lines(k))
    enddo
    if (mode.ge.0) then
      stop
    else
      call xit("regrid",mode)
    endif
  end subroutine usage

  !====================================================================================
  !--- Given src_grid and dst_grid containing grid description file names,
  !--- add definitions for the source grid and for the destination grid
  !--- This will also allocate and assign the associated mask and area arrays
  !---    src_grid_mask, dst_grid_mask, src_grid_area, dst_grid_area
  !====================================================================================
  subroutine define_grids()

    if ( len_trim(src_grid) < 1 .or. len_trim(dst_grid) < 1 ) then
      write(6,*)"define_grids: Both src_grid and dst_grid must be defined."
      call xit("define_grids",-1)
    endif

    !---------------------------------------------------------------------------------
    !--- Source Grid
    !---------------------------------------------------------------------------------

    !--- Add a definition for the source grid (internal to ccc_esmf)
    !--- This will define id_src_grid as well as
    !--- src_grid_dims, src_grid_edges, src_grid_size
    call add_grid(id_src_grid, src_grid, &
                  grid_dims=src_grid_dims, &
                  grid_edges=src_grid_edges, &
                  grid_size=src_grid_size)

    !--- Retrieve a pointer to the src grid mask
    !--- Note that reassigning src_grid_mask here will change the mask that is used in
    !--- the remap calculation below, since src_grid_mask is a pointer to that mask
    if ( associated(src_grid_mask) ) deallocate(src_grid_mask)
    call ESMF_GridGetItem(grid_list(id_src_grid), &
                          staggerloc=ESMF_STAGGERLOC_CENTER, &
                          itemflag=ESMF_GRIDITEM_MASK, &
                          farrayPtr=src_grid_mask, rc=rc)
    if (rc /= ESMF_SUCCESS) then
      write(6,*)'define_grids: Unable to retrive a pointer to the source grid mask.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif
    if ( .not. associated(src_grid_mask) ) then
      write(6,*)'define_grids: Pointer to the source grid mask is not associated.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif

    !--- Retrieve a pointer to the src grid area
    !--- Note that reassigning src_grid_area here will change the area internal to ccc_esmf
    if ( associated(src_grid_area) ) deallocate(src_grid_area)
    call ESMF_GridGetItem(grid_list(id_src_grid), &
                          staggerloc=ESMF_STAGGERLOC_CENTER, &
                          itemflag=ESMF_GRIDITEM_AREA, &
                          farrayPtr=src_grid_area, rc=rc)
    if (rc /= ESMF_SUCCESS) then
      write(6,*)'define_grids: Unable to retrive a pointer to the source grid area.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif
    if ( .not. associated(src_grid_area) ) then
      write(6,*)'define_grids: Pointer to the source grid area is not associated.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif

    !--- Retrieve a pointer to the src cell center lon (dim 1) coordinate values
    if ( associated(src_cent_lon) ) deallocate(src_cent_lon)
    call ESMF_GridGetCoord(grid_list(id_src_grid), coordDim=one, localDE=zero, &
                           staggerloc=ESMF_STAGGERLOC_CENTER, &
                           farrayPtr=src_cent_lon, rc=rc)
    if (rc /= ESMF_SUCCESS) then
      write(6,*)'define_grids: Unable to retrive a pointer to the source grid center lon.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif
    if ( .not. associated(src_cent_lon) ) then
      write(6,*)'define_grids: Pointer to the source grid center lon is not associated.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif

    !--- Retrieve a pointer to the src cell center lat (dim 2) coordinate values
    if ( associated(src_cent_lat) ) deallocate(src_cent_lat)
    call ESMF_GridGetCoord(grid_list(id_src_grid), coordDim=two, localDE=zero, &
                           staggerloc=ESMF_STAGGERLOC_CENTER, &
                           farrayPtr=src_cent_lat, rc=rc)
    if (rc /= ESMF_SUCCESS) then
      write(6,*)'define_grids: Unable to retrive a pointer to the source grid center lat.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif
    if ( .not. associated(src_cent_lat) ) then
      write(6,*)'define_grids: Pointer to the source grid center lat is not associated.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif

    !--- Retrieve a pointer to the src cell corner lon (dim 1) coordinate values
    if ( associated(src_corn_lon) ) deallocate(src_corn_lon)
    call ESMF_GridGetCoord(grid_list(id_src_grid), coordDim=one, localDE=zero, &
                           staggerloc=ESMF_STAGGERLOC_CORNER, &
                           farrayPtr=src_corn_lon, rc=rc)
    if (rc /= ESMF_SUCCESS) then
      write(6,*)'define_grids: Unable to retrive a pointer to the source grid corner lon.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif
    if ( .not. associated(src_corn_lon) ) then
      write(6,*)'define_grids: Pointer to the source grid corner lon is not associated.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif

    !--- Retrieve a pointer to the src cell corner lat (dim 2) coordinate values
    if ( associated(src_corn_lat) ) deallocate(src_corn_lat)
    call ESMF_GridGetCoord(grid_list(id_src_grid), coordDim=two, localDE=zero, &
                           staggerloc=ESMF_STAGGERLOC_CORNER, &
                           farrayPtr=src_corn_lat, rc=rc)
    if (rc /= ESMF_SUCCESS) then
      write(6,*)'define_grids: Unable to retrive a pointer to the source grid corner lat.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif
    if ( .not. associated(src_corn_lat) ) then
      write(6,*)'define_grids: Pointer to the source grid corner lat is not associated.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif

    write(6,*)'define_grids: Added grid ',trim(src_grid)

    !---------------------------------------------------------------------------------
    !--- Destination Grid
    !---------------------------------------------------------------------------------

    !--- Add a definition for the destination grid (internal to ccc_esmf)
    !--- This will define id_dst_grid as well as
    !--- dst_grid_dims, dst_grid_edges, dst_grid_size
    call add_grid(id_dst_grid, dst_grid, &
                  grid_dims=dst_grid_dims, &
                  grid_edges=dst_grid_edges, &
                  grid_size=dst_grid_size)

    !--- Retrieve a pointer to the dst grid mask
    !--- Note that reassigning dst_grid_mask here will change the mask that is used in
    !--- the remap calculation below, since dst_grid_mask is a pointer to that mask
    call ESMF_GridGetItem(grid_list(id_dst_grid), &
                          staggerloc=ESMF_STAGGERLOC_CENTER, &
                          itemflag=ESMF_GRIDITEM_MASK, &
                          farrayPtr=dst_grid_mask, rc=rc)
    if (rc /= ESMF_SUCCESS) then
      write(6,*)'define_grids: Unable to retrive a pointer to the source grid mask.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif
    if ( .not. associated(dst_grid_mask) ) then
      write(6,*)'define_grids: Pointer to the source grid mask is not associated.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif
    if ( size(dst_grid_mask,1) /= dst_grid_dims(1) ) then
      write(6,*)"define_grids: First dim of dst grid does not match dst_grid_dims."
      write(6,*)"   mask_dim(1),dst_grid_dims(1) = ",size(dst_grid_mask,1),dst_grid_dims(1)
      call xit("define_grids",-1)
    endif
    if ( size(dst_grid_mask,2) /= dst_grid_dims(2) ) then
      write(6,*)"define_grids: Second dim of dst grid does not match dst_grid_dims."
      write(6,*)"   mask_dim(2),dst_grid_dims(2) = ",size(dst_grid_mask,2),dst_grid_dims(2)
      call xit("define_grids",-1)
    endif

    !--- Retrieve a pointer to the dst grid area
    !--- Note that reassigning dst_grid_area here will change the area internal to ccc_esmf
    call ESMF_GridGetItem(grid_list(id_dst_grid), &
                          staggerloc=ESMF_STAGGERLOC_CENTER, &
                          itemflag=ESMF_GRIDITEM_AREA, &
                          farrayPtr=dst_grid_area, rc=rc)
    if (rc /= ESMF_SUCCESS) then
      write(6,*)'define_grids: Unable to retrive a pointer to the destination grid area.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif
    if ( .not. associated(dst_grid_area) ) then
      write(6,*)'define_grids: Pointer to the destination grid area is not associated.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif

    !--- Retrieve a pointer to the dst cell center lon (dim 1) coordinate values
    if ( associated(dst_cent_lon) ) deallocate(dst_cent_lon)
    call ESMF_GridGetCoord(grid_list(id_dst_grid), coordDim=one, localDE=zero, &
                           staggerloc=ESMF_STAGGERLOC_CENTER, &
                           farrayPtr=dst_cent_lon, rc=rc)
    if (rc /= ESMF_SUCCESS) then
      write(6,*)'define_grids: Unable to retrive a pointer to the destination grid center lon.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif
    if ( .not. associated(dst_cent_lon) ) then
      write(6,*)'define_grids: Pointer to the destination grid center lon is not associated.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif

    !--- Retrieve a pointer to the dst cell center lat (dim 2) coordinate values
    if ( associated(dst_cent_lat) ) deallocate(dst_cent_lat)
    call ESMF_GridGetCoord(grid_list(id_dst_grid), coordDim=two, localDE=zero, &
                           staggerloc=ESMF_STAGGERLOC_CENTER, &
                           farrayPtr=dst_cent_lat, rc=rc)
    if (rc /= ESMF_SUCCESS) then
      write(6,*)'define_grids: Unable to retrive a pointer to the destination grid center lat.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif
    if ( .not. associated(dst_cent_lat) ) then
      write(6,*)'define_grids: Pointer to the destination grid center lat is not associated.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif

    !--- Retrieve a pointer to the dst cell corner lon (dim 1) coordinate values
    if ( associated(dst_corn_lon) ) deallocate(dst_corn_lon)
    call ESMF_GridGetCoord(grid_list(id_dst_grid), coordDim=one, localDE=zero, &
                           staggerloc=ESMF_STAGGERLOC_CORNER, &
                           farrayPtr=dst_corn_lon, rc=rc)
    if (rc /= ESMF_SUCCESS) then
      write(6,*)'define_grids: Unable to retrive a pointer to the destination grid corner lon.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif
    if ( .not. associated(dst_corn_lon) ) then
      write(6,*)'define_grids: Pointer to the destination grid corner lon is not associated.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif

    !--- Retrieve a pointer to the dst cell corner lat (dim 2) coordinate values
    if ( associated(dst_corn_lat) ) deallocate(dst_corn_lat)
    call ESMF_GridGetCoord(grid_list(id_dst_grid), coordDim=two, localDE=zero, &
                           staggerloc=ESMF_STAGGERLOC_CORNER, &
                           farrayPtr=dst_corn_lat, rc=rc)
    if (rc /= ESMF_SUCCESS) then
      write(6,*)'define_grids: Unable to retrive a pointer to the destination grid corner lat.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif
    if ( .not. associated(dst_corn_lat) ) then
      write(6,*)'define_grids: Pointer to the destination grid corner lat is not associated.'
      call ESMF_Finalize(endflag=ESMF_END_ABORT)
      call xit("define_grids",-2)
    endif

    write(6,*)'define_grids: Added grid ',trim(dst_grid)

  end subroutine define_grids

  !====================================================================================
  !--- Create fractional land mask on AGCM grid from NEMO land mask
  !====================================================================================
  subroutine create_fractional_land_mask()

    !--- Local
    integer :: idx1, ix, iy
    integer :: regridID1, regridID2
    real(ESMF_KIND_R8), pointer :: cell_frac(:,:)
    character(256) :: out_file = " "

    !--- Define source and destination grids
    call define_grids()

    !--- Determine which grid masks are to be used
    if ( associated(src_mask_values) ) deallocate(src_mask_values)
    allocate( src_mask_values(1) )
    src_mask_values = 0
    if ( associated(dst_mask_values) ) deallocate(dst_mask_values)
    allocate( dst_mask_values(1) )
    dst_mask_values = 0

    !--- Set regrid method and normalization option here
    regrid_method = "conservative"
    this_norm_opt = "fracarea"

    !--- This invocation of find_regrid_id will cause weights to be generated
    !--- and a new regrid operation defined in regrid_list(regridID1)
    use_src_mask=.false.
    use_dst_mask=.true.
    call find_regrid_id(regridID1, src_grid, dst_grid, regrid_method,  &
                        use_src_mask, use_dst_mask, &
                        src_mask_values, dst_mask_values, &
                        norm_type=this_norm_opt, create=.true.)

    !--- This call to find_regrid_id will assign the fraction of each source
    !--- and destination grid cell that participates in the regridding to
    !---    regrid_list(regridID1)%src_frac(:,:)
    !---    regrid_list(regridID1)%dst_frac(:,:)
    if ( .not. associated(regrid_list(regridID1)%src_frac) ) then
      write(6,*)"create_fractional_land_mask: src_frac is not defined."
      call xit("create_fractional_land_mask",-1)
    endif
    if ( .not. associated(regrid_list(regridID1)%dst_frac) ) then
      write(6,*)"create_fractional_land_mask: dst_frac is not defined."
      call xit("create_fractional_land_mask",-1)
    endif

    !--- This invocation of find_regrid_id will cause weights to be generated
    !--- and a new regrid operation defined in regrid_list(regridID2)
    use_src_mask=.true.
    use_dst_mask=.false.
    call find_regrid_id(regridID2, dst_grid, src_grid, regrid_method,  &
                        use_src_mask, use_dst_mask, &
                        src_mask_values, dst_mask_values, &
                        norm_type=this_norm_opt, create=.true.)

    !--- This call to find_regrid_id will assign the fraction of each source
    !--- and destination grid cell that participates in the regridding to
    !---    regrid_list(regridID2)%src_frac(:,:)
    !---    regrid_list(regridID2)%dst_frac(:,:)
    if ( .not. associated(regrid_list(regridID2)%src_frac) ) then
      write(6,*)"create_fractional_land_mask: src_frac is not defined."
      call xit("create_fractional_land_mask",-1)
    endif
    if ( .not. associated(regrid_list(regridID2)%dst_frac) ) then
      write(6,*)"create_fractional_land_mask: dst_frac is not defined."
      call xit("create_fractional_land_mask",-1)
    endif

    !--- Print coordinate values where src_frac > 1
    !--- This happens along the north fold where src_frac could be as much a 2.0
    do ix=1,size(regrid_list(regridID1)%src_frac, dim=1)
      do iy=1,size(regrid_list(regridID1)%src_frac, dim=2)
        if ( regrid_list(regridID1)%src_frac(ix,iy) > 1.00000000001 ) then
          write(6,*)"ix,iy: ",ix,iy, &
                    " lon,lat: ",src_cent_lon(ix,iy),src_cent_lat(ix,iy), &
                    " src_frac: ",regrid_list(regridID1)%src_frac(ix,iy)
        endif
      enddo
    enddo

    !--- Write the src fraction stored in regrid_list(regridID1)%src_frac(:,:) to disk
    open(701, file="frac_land_mask_file_all", form='unformatted')
      ibuf(1:8) = 1
      ibuf(1) = transfer("GRID", 1)
      ibuf(3) = transfer("SRCF", 1)
      ibuf(5) = src_grid_dims(1)
      ibuf(6) = src_grid_dims(2)
      ibuf(7) = 0

      !--- The regrid operation T63 to ORCA1
      wrk(1:ibuf(5)*ibuf(6)) = reshape(regrid_list(regridID1)%src_frac, (/ibuf(5)*ibuf(6)/))
      !--- In the first version of this fractional land mask data, all cell fractions
      !--- that were greater than 1 were set to 1
      where( wrk(1:ibuf(5)*ibuf(6)) > 1.0 ) wrk(1:ibuf(5)*ibuf(6)) = 1.0
      call putfld2(701,wrk,ibuf,maxx)

      !--- The inverse regrid operation ORCA1 to T63
      ibuf(3) = transfer("DSTF", 1)
      wrk(1:ibuf(5)*ibuf(6)) = reshape(regrid_list(regridID2)%dst_frac, (/ibuf(5)*ibuf(6)/))
      !--- Set all cell fractions > 1 to 1
      where( wrk(1:ibuf(5)*ibuf(6)) > 1.0 ) wrk(1:ibuf(5)*ibuf(6)) = 1.0
      call putfld2(701,wrk,ibuf,maxx)

      !--- DBG output
      ibuf(3) = transfer("OBPT", 1)
      wrk(1:ibuf(5)*ibuf(6)) = reshape(regrid_list(regridID1)%src_frac, (/ibuf(5)*ibuf(6)/))
      !--- Set all cell fractions <= 1 to 0
      where( wrk(1:ibuf(5)*ibuf(6)) <= 1.0 ) wrk(1:ibuf(5)*ibuf(6)) = 0.0
      call putfld2(701,wrk,ibuf,maxx)

      !--- The difference in T63 cell fraction between T63 to ORCA1 and ORCA1 to T63
      ibuf(3) = transfer("DIFF", 1)
      wrk(1:ibuf(5)*ibuf(6)) = reshape(regrid_list(regridID1)%src_frac, (/ibuf(5)*ibuf(6)/)) &
                              -reshape(regrid_list(regridID2)%dst_frac, (/ibuf(5)*ibuf(6)/))
      call putfld2(701,wrk,ibuf,maxx)

      write(6,'(a,a4,1x,i10,1x,a4,1x,5i8)')"Created frac_land_mask_file_all"
    close(701)

    !--- Write another file containing only the T63 cell fraction
    ibuf(1:8) = 1
    ibuf(1) = transfer("GRID", 1)
    ibuf(3) = transfer("SRCF", 1)
    ibuf(7) = 0

    !--- The regrid operation T63 to ORCA1
    if ( frac_land_mask_add_olap ) then
      !--- Include the overlap longitude in the output file
      out_file = "t63_to_orca1_grid_cell_fraction_129x64"
      open(701, file=trim(out_file), form='unformatted')
      ibuf(5) = src_grid_dims(1) + 1
      ibuf(6) = src_grid_dims(2)
      if ( associated(cell_frac) ) deallocate(cell_frac)
      allocate( cell_frac(ibuf(5),ibuf(6)) )
      cell_frac(1:ibuf(5)-1,1:ibuf(6)) = regrid_list(regridID1)%src_frac(1:ibuf(5)-1,1:ibuf(6))
      cell_frac(ibuf(5),:) = cell_frac(1,:)
    else
      !--- Do not include the overlap longitude in the output file
      out_file = "t63_to_orca1_grid_cell_fraction_128x64"
      open(701, file=trim(out_file), form='unformatted')
      ibuf(5) = src_grid_dims(1)
      ibuf(6) = src_grid_dims(2)
      if ( associated(cell_frac) ) deallocate(cell_frac)
      allocate( cell_frac(ibuf(5),ibuf(6)) )
      cell_frac(1:ibuf(5),1:ibuf(6)) = regrid_list(regridID1)%src_frac(1:ibuf(5),1:ibuf(6))
    endif
    wrk(1:ibuf(5)*ibuf(6)) = reshape(cell_frac, (/ibuf(5)*ibuf(6)/))
    deallocate(cell_frac)
    !--- In the first version of this fractional land mask data, all cell fractions
    !--- that were greater than 1 were set to 1
    where( wrk(1:ibuf(5)*ibuf(6)) > 1.0 ) wrk(1:ibuf(5)*ibuf(6)) = 1.0
    call putfld2(701,wrk,ibuf,maxx)
    write(6,'(3a,a4,1x,i10,1x,a4,1x,5i8)')"Created ",trim(out_file)," containing ",ibuf(1:8)
    close(701)

    if ( .false. ) then

    !--- Define cpl_vinfo_t variables (internal to ccc_esmf)
    allocate(src_vinfo)
    src_vinfo%name     = "src_field"
    src_vinfo%grid     = trim(src_grid)
    src_vinfo%tag      = cpl_var_new_tag(5000)
    src_vinfo%size     = src_grid_dims(1)*src_grid_dims(2)
    src_vinfo%nlon     = src_grid_dims(1)  !--- Do not include overlap
    src_vinfo%nlat     = src_grid_dims(2)
    src_vinfo%olap     = 0
    src_vinfo%ncopy    = 1
    src_vinfo%freq     = 0
    src_vinfo%ccc_name = " SRC"
    src_vinfo%regrid%method   = trim(regrid_method)
    src_vinfo%regrid%norm_opt = trim(this_norm_opt)
    src_vinfo%regrid%masked   = use_src_mask
    !--- Add this variable defintion, overwriting any existing def with the same name
    call add_cpl_var(src_vinfo, exists=1)

    allocate(dst_vinfo)
    dst_vinfo%name     = "dst_field"
    dst_vinfo%grid     = trim(dst_grid)
    dst_vinfo%tag      = cpl_var_new_tag(5000)
    dst_vinfo%size     = dst_grid_dims(1)*dst_grid_dims(2)
    dst_vinfo%nlon     = dst_grid_dims(1)  !--- Do not include overlap
    dst_vinfo%nlat     = dst_grid_dims(2)
    dst_vinfo%olap     = 0
    dst_vinfo%ncopy    = 1
    dst_vinfo%freq     = 0
    dst_vinfo%ccc_name = " DST"
    dst_vinfo%regrid%method   = trim(regrid_method)
    dst_vinfo%regrid%norm_opt = trim(this_norm_opt)
    dst_vinfo%regrid%masked   = use_dst_mask
    !--- Add this variable defintion, overwriting any existing def with the same name
    call add_cpl_var(dst_vinfo, exists=1)

    if ( associated(src_data)  ) deallocate(src_data)
    if ( associated(dst_data)  ) deallocate(dst_data)
    if ( associated(itmp_data) ) deallocate(itmp_data)

    allocate( src_data(src_grid_dims(1)*src_grid_dims(2)) )
    allocate( itmp_data(src_grid_dims(1)*src_grid_dims(2)) )
    src_data = 2.0
    itmp_data = reshape( src_grid_mask, (/src_grid_dims(1)*src_grid_dims(2)/) )
    where ( itmp_data == 0 ) src_data = 1.0

    !--- Make this a bilinear regrid operation
    dst_vinfo%regrid%method   = "bilinear"
    call add_cpl_var(dst_vinfo, exists=1)

    call regrid_Field(src_vinfo%name, src_data, dst_vinfo%name, dst_data, &
                      use_src_mask=.false., use_dst_mask=.true. )

    !--- Reset the regrid method in dst_vinfo
    dst_vinfo%regrid%method   = trim(regrid_method)
    call add_cpl_var(dst_vinfo, exists=1)

    ibuf(1:8) = 1
    ibuf(1) = transfer("GRID", 1)
    ibuf(7) = 0
    open(701, file="frac_land_mask_file", form='unformatted')
      ibuf(3) = transfer("FRAC", 1)
      ibuf(5) = dst_grid_dims(1)
      ibuf(6) = dst_grid_dims(2)
      wrk(1:ibuf(5)*ibuf(6)) = dst_data(1:ibuf(5)*ibuf(6))
      call putfld2(701,wrk,ibuf,maxx)
      write(6,'(a,a4,1x,i10,1x,a4,1x,5i8)')"Wrote to frac_land_mask_file: ",ibuf(1:8)
    close(701)

    if ( associated(src_data)  ) deallocate(src_data)
    if ( associated(dst_data)  ) deallocate(dst_data)
    if ( associated(itmp_data) ) deallocate(itmp_data)

    endif

  end subroutine create_fractional_land_mask

  !====================================================================================
  !--- testing
  !====================================================================================
  subroutine test_grid()
    !--- Local
    real(kind=8), pointer :: src_data2d(:,:)
    real(kind=8), pointer :: dst_data2d(:,:)
    real(kind=8), pointer :: src_lon(:), src_lat(:), dst_lon(:), dst_lat(:)
    real(kind=8), pointer :: src_lon_bnds(:), src_lat_bnds(:), dst_lon_bnds(:), dst_lat_bnds(:)
    integer :: nx_src, ny_src, nx_dst, ny_dst, ix, iy, idx, member, k_test_func
    integer :: nx_src_frac, ny_src_frac, nx_dst_frac, ny_dst_frac
    logical :: use_src_frac, use_dst_frac, renorm_destarea
    integer :: curr_src_grid_id, curr_dst_grid_id, err2_regrid_id
    real(ESMF_KIND_R8),    pointer :: err1_data(:,:), err2_data(:,:)
    real(ESMF_KIND_R8),    pointer :: src_area(:,:), dst_area(:,:)
    integer(ESMF_KIND_I4), pointer :: src_mask(:,:), dst_mask(:,:)
    integer, pointer :: curr_src_mask(:,:), curr_dst_mask(:,:)
    real(ESMF_KIND_R8),    pointer :: src_frac(:,:), dst_frac(:,:)
    real(ESMF_KIND_R8) :: src_avg, dst_avg
    integer(ESMF_KIND_I), parameter ::  one=1
    integer(ESMF_KIND_I), parameter ::  two=2

    integer(kind=4) :: fid, rec_dim_id, rec_id
    integer(kind=4) :: xsrc_id, ysrc_id, xsrc_dim_id, ysrc_dim_id
    integer(kind=4) :: xdst_id, ydst_id, xdst_dim_id, ydst_dim_id
    integer(kind=4) :: xsrc_bnds_id, ysrc_bnds_id, xsrc_bnds_dim_id, ysrc_bnds_dim_id
    integer(kind=4) :: xdst_bnds_id, ydst_bnds_id, xdst_bnds_dim_id, ydst_bnds_dim_id
    integer(kind=4) :: src_data_id, dst_data_id, err1_data_id, err2_data_id
    integer(kind=4) :: src_mask_id, dst_mask_id
    integer(kind=4) :: src_area_id, dst_area_id
    integer(kind=4) :: src_frac_id, dst_frac_id

    integer(kind=4) :: descll_id, descal_id, description_id
    character(80) :: description(8)

    !--- The time associated with this variable
    real(kind=8) :: rec_val

    !--- The name of variable to be written
    character(32) :: var_name

    !--- The name of the output file
    character(128) :: file_out_name, strng1, strng2, strng3

    !--- Add a definition for the source grid (internal to ccc_esmf)
    call add_grid(id_src_grid, src_grid, &
                  grid_dims=src_grid_dims, &
                  grid_edges=src_grid_edges, &
                  grid_size=src_grid_size, &
                  cartesian=.true.)
    nx_src = src_grid_dims(1)
    ny_src = src_grid_dims(2)

    !--- Add a definition for the destination grid (internal to ccc_esmf)
    call add_grid(id_dst_grid, dst_grid, &
                  grid_dims=dst_grid_dims, &
                  grid_edges=dst_grid_edges, &
                  grid_size=dst_grid_size, &
                  cartesian=.true.)
    nx_dst = dst_grid_dims(1)
    ny_dst = dst_grid_dims(2)

    !--- Get the x and y center coordinate values on the src and dst grid
    call ESMF_GridGetCoord(grid_list(id_src_grid), coordDim=one, &
                           staggerloc=ESMF_STAGGERLOC_CENTER, &
                           farrayPtr=src_lon, rc=rc)
    lnum = __LINE__ ; call handle_ESMF_err(rc, &
        "Cannot retrieve src grid X values", line=lnum, method=method)

    call ESMF_GridGetCoord(grid_list(id_src_grid), coordDim=two, &
                           staggerloc=ESMF_STAGGERLOC_CENTER, &
                           farrayPtr=src_lat, rc=rc)
    lnum = __LINE__ ; call handle_ESMF_err(rc, &
        "Cannot retrieve src grid Y values", line=lnum, method=method)

    call ESMF_GridGetCoord(grid_list(id_dst_grid), coordDim=one, &
                           staggerloc=ESMF_STAGGERLOC_CENTER, &
                           farrayPtr=dst_lon, rc=rc)
    lnum = __LINE__ ; call handle_ESMF_err(rc, &
        "Cannot retrieve dst grid X values", line=lnum, method=method)

    call ESMF_GridGetCoord(grid_list(id_dst_grid), coordDim=two, &
                           staggerloc=ESMF_STAGGERLOC_CENTER, &
                           farrayPtr=dst_lat, rc=rc)
    lnum = __LINE__ ; call handle_ESMF_err(rc, &
        "Cannot retrieve dst grid Y values", line=lnum, method=method)

    if ( verbose > 1 ) then
      write(6,*)"src_lon: ",src_lon
      write(6,*)"src_lat: ",src_lat
      write(6,*)"dst_lon: ",dst_lon
      write(6,*)"dst_lat: ",dst_lat
    endif

    !--- Get the x and y corner coordinate values on the src and dst grid
    call ESMF_GridGetCoord(grid_list(id_src_grid), coordDim=one, &
                           staggerloc=ESMF_STAGGERLOC_CORNER, &
                           farrayPtr=src_lon_bnds, rc=rc)
    lnum = __LINE__ ; call handle_ESMF_err(rc, &
        "Cannot retrieve src grid X bounds values", line=lnum, method=method)

    call ESMF_GridGetCoord(grid_list(id_src_grid), coordDim=two, &
                           staggerloc=ESMF_STAGGERLOC_CORNER, &
                           farrayPtr=src_lat_bnds, rc=rc)
    lnum = __LINE__ ; call handle_ESMF_err(rc, &
        "Cannot retrieve src grid Y bounds values", line=lnum, method=method)

    call ESMF_GridGetCoord(grid_list(id_dst_grid), coordDim=one, &
                           staggerloc=ESMF_STAGGERLOC_CORNER, &
                           farrayPtr=dst_lon_bnds, rc=rc)
    lnum = __LINE__ ; call handle_ESMF_err(rc, &
        "Cannot retrieve dst grid X bounds values", line=lnum, method=method)

    call ESMF_GridGetCoord(grid_list(id_dst_grid), coordDim=two, &
                           staggerloc=ESMF_STAGGERLOC_CORNER, &
                           farrayPtr=dst_lat_bnds, rc=rc)
    lnum = __LINE__ ; call handle_ESMF_err(rc, &
        "Cannot retrieve dst grid Y bounds values", line=lnum, method=method)

    if ( verbose > 1 ) then
      write(6,*)"src_lon_bnds: ",src_lon_bnds
      write(6,*)"src_lat_bnds: ",src_lat_bnds
      write(6,*)"dst_lon_bnds: ",dst_lon_bnds
      write(6,*)"dst_lat_bnds: ",dst_lat_bnds
    endif

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

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

    !--- Retrieve a pointer to the grid mask in the src grid
    call ESMF_GridGetItem(grid_list(id_src_grid), &
                          staggerloc=ESMF_STAGGERLOC_CENTER, &
                          itemflag=ESMF_GRIDITEM_MASK, &
                          farrayPtr=src_mask, rc=rc)
    lnum = __LINE__ ; call handle_ESMF_err(rc, &
        "Cannot retrieve src grid mask", line=lnum, method=method)

    !--- Retrieve a pointer to the grid mask in the dst grid
    call ESMF_GridGetItem(grid_list(id_dst_grid), &
                          staggerloc=ESMF_STAGGERLOC_CENTER, &
                          itemflag=ESMF_GRIDITEM_MASK, &
                          farrayPtr=dst_mask, rc=rc)
    lnum = __LINE__ ; call handle_ESMF_err(rc, &
        "Cannot retrieve dst grid mask", line=lnum, method=method)

    !--- Determine which grid masks are to be used
    use_src_mask=.true.
    use_dst_mask=.true.

    !--- Determine the value in each mask that indicates that that cell
    !--- is not to be used in the remapping operation
    !--- As defined here, a cell mask value = 0 indicates that cell is not used
    if ( associated(src_mask_values) ) deallocate(src_mask_values)
    allocate( src_mask_values(1) )
    src_mask_values = 0
    if ( associated(dst_mask_values) ) deallocate(dst_mask_values)
    allocate( dst_mask_values(1) )
    dst_mask_values = 0

    write(6,*)"regrid: src_mask_values = ",src_mask_values
    write(6,*)"regrid: dst_mask_values = ",dst_mask_values

    !--- Allocate space for the source and destination data
    if ( associated(src_data)  ) deallocate(src_data)
    if ( associated(dst_data)  ) deallocate(dst_data)
    if ( associated(src_data2d)  ) deallocate(src_data2d)
    if ( associated(dst_data2d)  ) deallocate(dst_data2d)
    allocate( src_data(nx_src*ny_src) )
    allocate( dst_data(nx_dst*ny_dst) )
    allocate( src_data2d(nx_src,ny_src) )
    allocate( dst_data2d(nx_dst,ny_dst) )

    !--- Define cpl_vinfo_t variables (internal to ccc_esmf)
    allocate(src_vinfo)
    src_vinfo%name     = "src_field"
    src_vinfo%grid     = trim(src_grid)
    src_vinfo%tag      = cpl_var_new_tag(5000)
    src_vinfo%size     = src_grid_dims(1)*src_grid_dims(2)
    src_vinfo%nlon     = src_grid_dims(1)  !--- Do not include overlap
    src_vinfo%nlat     = src_grid_dims(2)
    src_vinfo%olap     = 0
    src_vinfo%ncopy    = 1
    src_vinfo%freq     = 0
    src_vinfo%ccc_name = " SRC"
    src_vinfo%regrid%method   = trim(regrid_method)
    src_vinfo%regrid%norm_opt = trim(this_norm_opt)
    src_vinfo%regrid%masked   = use_src_mask
    !--- Add this variable defintion, overwriting any existing def with the same name
    call add_cpl_var(src_vinfo, exists=1)

    allocate(dst_vinfo)
    dst_vinfo%name     = "dst_field"
    dst_vinfo%grid     = trim(dst_grid)
    dst_vinfo%tag      = cpl_var_new_tag(5000)
    dst_vinfo%size     = dst_grid_dims(1)*dst_grid_dims(2)
    dst_vinfo%nlon     = dst_grid_dims(1)  !--- Do not include overlap
    dst_vinfo%nlat     = dst_grid_dims(2)
    dst_vinfo%olap     = 0
    dst_vinfo%ncopy    = 1
    dst_vinfo%freq     = 0
    dst_vinfo%ccc_name = " DST"
    dst_vinfo%regrid%method   = trim(regrid_method)
    dst_vinfo%regrid%norm_opt = trim(this_norm_opt)
    dst_vinfo%regrid%masked   = use_dst_mask
    !--- Add this variable defintion, overwriting any existing def with the same name
    call add_cpl_var(dst_vinfo, exists=1)

    !--- Assign source data values
    if ( nx_src == 3 .and. ny_src == 3 ) then
      k_test_func = 5
    else if ( nx_src == 5 .and. ny_src == 5 ) then
      k_test_func = 6
    else if ( nx_src == 6 .and. ny_src == 6 .and. interp_mask ) then
      k_test_func = 7
    else if ( nx_src == 6 .and. ny_src == 6 .and. interp_sicn ) then
      k_test_func = 8
    else
      k_test_func = 2
    endif
    do ix=1,nx_src
      do iy=1,ny_src
        src_data2d(ix,iy) = test_function( k_test_func, src_lon(ix), src_lat(iy) )
      enddo
    enddo
    src_data = reshape(src_data2d, (/nx_src*ny_src/))

    !--- Define the output file name
    strng1 = " "
    if ( src_grid(1:10) .eq. "grid_desc_" ) then
      strng1 = trim(src_grid(11:))
    else
      strng1 = trim(src_grid)
    endif
    if ( strng1(1:4) .eq. "test" ) then
      strng1 = trim(strng1(5:))
    endif
    ix = len_trim(strng1)
    if ( strng1(ix-2:ix) .eq. ".nc" ) then
      strng1 = strng1(1:ix-3)
    endif

    strng2 = " "
    if ( dst_grid(1:10) .eq. "grid_desc_" ) then
      strng2 = trim(dst_grid(11:))
    else
      strng2 = trim(dst_grid)
    endif
    if ( strng2(1:4) .eq. "test" ) then
      strng2 = trim(strng2(5:))
    endif
    ix = len_trim(strng2)
    if ( strng2(ix-2:ix) .eq. ".nc" ) then
      strng2 = strng2(1:ix-3)
    endif

    strng3 = " "
    write(strng3,'(a,i2.2)')"func",k_test_func

    file_out_name = " "
    file_out_name = "regrid_"//trim(strng3)//"_"//trim(strng1)//"_to_"//trim(strng2)//".nc"

    !--- Create a netcdf file to contain the test data
    call ncdf_open(fid, file_out_name)

    !--- Add coordinate variables
    call ncdf_add_coord(fid, xsrc_id, xsrc_dim_id, "xsrc", vals=src_lon)
    call ncdf_add_coord(fid, ysrc_id, ysrc_dim_id, "ysrc", vals=src_lat)
    call ncdf_add_coord(fid, xdst_id, xdst_dim_id, "xdst", vals=dst_lon)
    call ncdf_add_coord(fid, ydst_id, ydst_dim_id, "ydst", vals=dst_lat)

    !--- Add coordinate bounds
    call ncdf_add_coord(fid, xsrc_bnds_id, xsrc_bnds_dim_id, "xsrc_bnds", vals=src_lon_bnds)
    call ncdf_add_coord(fid, ysrc_bnds_id, ysrc_bnds_dim_id, "ysrc_bnds", vals=src_lat_bnds)
    call ncdf_add_coord(fid, xdst_bnds_id, xdst_bnds_dim_id, "xdst_bnds", vals=dst_lon_bnds)
    call ncdf_add_coord(fid, ydst_bnds_id, ydst_bnds_dim_id, "ydst_bnds", vals=dst_lat_bnds)

    !--- Define a record dimension named "member"
    !--- NOTE: Apparently the record dimension cannot be the first thing defined in the file
    call ncdf_add_dim(fid, rec_id, "member", record_dim=.true.)

    !--- Dimensions not associated with a variable must be
    !--- defined after the record dimension is defined

    !--- Add other variable definitions
    call ncdf_add_var(fid, src_data_id, "src_data", (/xsrc_id,ysrc_id,rec_id/))
    call ncdf_add_att(fid, src_data_id, "long_name", "Source data")

    call ncdf_add_var(fid, dst_data_id, "dst_data", (/xdst_id,ydst_id,rec_id/))
    call ncdf_add_att(fid, dst_data_id, "long_name", "Interpolated data")

    call ncdf_add_var(fid, err1_data_id, "err1_data", (/xdst_id,ydst_id,rec_id/))
    call ncdf_add_att(fid, err1_data_id, "long_name", "Error in interpolated data")

    call ncdf_add_var(fid, err2_data_id, "err2_data", (/xsrc_id,ysrc_id,rec_id/))
    call ncdf_add_att(fid, err2_data_id, "long_name", "Error in re-interpolated data")

    call ncdf_add_var(fid, src_mask_id, "src_mask", (/xsrc_id,ysrc_id,rec_id/), xtype=NF90_int)
    call ncdf_add_att(fid, src_mask_id, "long_name", "Source mask")
    call ncdf_write_var(fid, src_mask_id, vals=src_mask)

    call ncdf_add_var(fid, dst_mask_id, "dst_mask", (/xdst_id,ydst_id,rec_id/), xtype=NF90_int)
    call ncdf_add_att(fid, dst_mask_id, "long_name", "Destination mask")
    call ncdf_write_var(fid, dst_mask_id, vals=dst_mask)

    call ncdf_add_var(fid, src_frac_id, "src_frac", (/xsrc_id,ysrc_id,rec_id/))
    call ncdf_add_att(fid, src_frac_id, "long_name", "Source fraction")
    call ncdf_add_att(fid, src_frac_id, "comment", &
         "The fraction of each source cell participating in the remap")
    call ncdf_add_att(fid, src_frac_id, "caveat", &
         "src_frac is only defined for a conservative remap")

    call ncdf_add_var(fid, dst_frac_id, "dst_frac", (/xdst_id,ydst_id,rec_id/))
    call ncdf_add_att(fid, dst_frac_id, "long_name", "Destination fraction")
    call ncdf_add_att(fid, dst_frac_id, "comment", &
         "The fraction of each destination cell participating in the remap")
    call ncdf_add_att(fid, dst_frac_id, "caveat", &
         "dst_frac is only defined for a conservative remap")

    call ncdf_add_var(fid, src_area_id, "src_area", (/xsrc_id,ysrc_id/), vals=src_area)
    call ncdf_add_att(fid, src_area_id, "long_name", "Source area")

    call ncdf_add_var(fid, dst_area_id, "dst_area", (/xdst_id,ydst_id/), vals=dst_area)
    call ncdf_add_att(fid, dst_area_id, "long_name", "Destination area")

    !--- This is an example of how to define and assign a 1D character array in a netcdf file
    !--- The first (fastest changing) dimension must always be the length of a single element
    !--- in the character array. The second dimension will be the length of the array.
    !--- If the variable is to contain the record dimension then it must be the
    !--- last (slowest varying) dimension

    !--- Define a character array to hold a description of each test
    call ncdf_add_dim(fid, descll_id, "desc_ll", len(description(1)))
    call ncdf_add_dim(fid, descal_id, "desc_al", size(description))
    call ncdf_add_var(fid, description_id, "description", (/descll_id,descal_id,rec_id/), xtype=NF90_char)
    call ncdf_add_att(fid, description_id, "long_name", "Test descriptions")

    !--- Write global meta data
    call ncdf_add_att(fid, NF90_GLOBAL, "title", "ESMF Interpolation Tests")

    select case (k_test_func)
      case (1)
        !--- A constant value of one
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function", &
             "Constant value of 1")

      case (2)
        !--- Quadratic in x and y, constant on boundaries
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function", &
             "1 + 0.5*(x0*(x-x1)*(x2-x) + y0*(y-y1)*(y2-y))")
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function_x1_x2", (/test_func_xmin,test_func_xmax/))
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function_y1_y2", (/test_func_ymin,test_func_ymax/))
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function_x0", "4/(x2-x1)**2")
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function_y0", "4/(y2-y1)**2")

      case (3)
        !--- Bimodal trig function
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function", &
             "1 + cos(x0*(x-x1)*(x2-x)*pi)*cos(y0*(y-y1)*(y2-y)*pi)")
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function_x1_x2", (/test_func_xmin,test_func_xmax/))
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function_y1_y2", (/test_func_ymin,test_func_ymax/))
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function_x0", "4/(x2-x1)**2")
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function_y0", "4/(y2-y1)**2")

      case (4)
        !--- Higher frequency trig function
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function", &
             "1 + cos(x0*(x-x1)*(x2-x)*5*pi)*cos(y0*(y-y1)*(y2-y)*5*pi)")
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function_x1_x2", (/test_func_xmin,test_func_xmax/))
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function_y1_y2", (/test_func_ymin,test_func_ymax/))
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function_x0", "4/(x2-x1)**2")
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function_y0", "4/(y2-y1)**2")

      case (5)
         !--- Checker board pattern on a 3x3 grid
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function", &
             "Checker board pattern (1/10) on a 3x3 grid")

      case (6)
         !--- Checker board pattern on a 5x5 grid
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function", &
             "Checker board pattern (1/10) on a 5x5 grid")

      case (7)
         !--- Example of a land Mask on a 6x6 grid
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function", &
             "Example of a land Mask on a 6x6 grid")

      case (8)
         !--- Example of SICN on a 6x6 grid
        call ncdf_add_att(fid, NF90_GLOBAL, "test_function", &
             "Example of SICN on a 6x6 grid")

    end select

    member=0

    if ( associated(curr_src_mask) ) deallocate(curr_src_mask)
    allocate( curr_src_mask(nx_src,ny_src) )
    curr_src_mask = src_mask
    curr_src_grid_id = id_src_grid
    if ( nx_src == 6 .and. ny_src == 6 .and. k_test_func /= 7 ) then
      curr_src_mask(1:6,6) = (/ 1, 1, 1, 1, 1, 1 /)
      curr_src_mask(1:6,5) = (/ 0, 0, 0, 0, 1, 1 /)
      curr_src_mask(1:6,4) = (/ 1, 0, 0, 1, 1, 1 /)
      curr_src_mask(1:6,3) = (/ 1, 1, 1, 1, 1, 1 /)
      curr_src_mask(1:6,2) = (/ 1, 1, 1, 1, 0, 0 /)
      curr_src_mask(1:6,1) = (/ 1, 1, 1, 1, 0, 0 /)
      call add_grid(curr_src_grid_id, src_grid, imask=curr_src_mask, cartesian=.true.)
    endif

    if ( associated(curr_dst_mask) ) deallocate(curr_dst_mask)
    allocate( curr_dst_mask(nx_dst,ny_dst) )
    curr_dst_mask = dst_mask
    curr_dst_grid_id = id_dst_grid
    if ( nx_dst == 6 .and. ny_dst == 6 ) then
      curr_dst_mask(1:6,6) = (/ 1, 1, 1, 1, 1, 1 /)
      curr_dst_mask(1:6,5) = (/ 0, 0, 0, 0, 1, 1 /)
      curr_dst_mask(1:6,4) = (/ 1, 0, 0, 1, 1, 1 /)
      curr_dst_mask(1:6,3) = (/ 1, 1, 1, 1, 1, 1 /)
      curr_dst_mask(1:6,2) = (/ 1, 1, 1, 1, 0, 0 /)
      curr_dst_mask(1:6,1) = (/ 1, 1, 1, 1, 0, 0 /)
      call add_grid(curr_dst_grid_id, dst_grid, imask=curr_dst_mask, cartesian=.true.)
    endif


    do idx=1,10

      !--- Make test specific changes
      description(:) = " "
      select case (idx)
        case (1)
          member = member + 1
          regrid_method = "conservative"
          this_norm_opt = "fracarea"
          renorm_destarea = .true.
          use_src_mask = .true.
          use_dst_mask = .true.

          description(1) = "method = conservative"
          description(2) = "Norm type = fracarea"
          strng = " "
          write(strng,'(a,l1)')"Renormalize destarea = ",renorm_destarea
          description(3) = trim(strng)

        case (2)
          member = member + 1
          regrid_method = "conservative"
          this_norm_opt = "destarea"
          renorm_destarea = .true.
          use_src_mask = .true.
          use_dst_mask = .true.

          description(1) = "method = conservative"
          description(2) = "Norm type = destarea"
          strng = " "
          write(strng,'(a,l1)')"Renormalize destarea = ",renorm_destarea
          description(3) = trim(strng)

        case (3)
          member = member + 1
          regrid_method = "conservative"
          this_norm_opt = "destarea"
          renorm_destarea = .false.
          use_src_mask = .true.
          use_dst_mask = .true.

          description(1) = "method = conservative"
          description(2) = "Norm type = destarea"
          strng = " "
          write(strng,'(a,l1)')"Renormalize destarea = ",renorm_destarea
          description(3) = trim(strng)

        case (104)
          member = member + 1
          regrid_method = "conservative"
          this_norm_opt = "fracarea"
          renorm_destarea = .true.
          use_src_mask = .true.
          use_dst_mask = .true.

          !--- Change the grid mask and add a new definition for the source grid
          if ( associated(curr_src_mask) ) deallocate(curr_src_mask)
          allocate( curr_src_mask(nx_src,ny_src) )
          curr_src_mask = src_mask
          curr_src_mask(1,1) = 0
          call add_grid(curr_src_grid_id, src_grid, imask=curr_src_mask, cartesian=.true.)

          description(1) = "method = conservative"
          description(2) = "Norm type = fracarea"
          strng = " "
          write(strng,'(a,l1)')"Renormalize destarea = ",renorm_destarea
          description(3) = trim(strng)

        case (105)
          member = member + 1
          regrid_method = "conservative"
          this_norm_opt = "fracarea"
          renorm_destarea = .true.
          use_src_mask = .true.
          use_dst_mask = .true.

          !--- Reset the src mask to its original value
          !--- curr_src_mask = src_mask
          !--- call add_grid(curr_src_grid_id, src_grid, imask=curr_src_mask, cartesian=.true.)

          !--- Change the grid mask and add a new definition for the destination grid
          if ( associated(curr_dst_mask) ) deallocate(curr_dst_mask)
          allocate( curr_dst_mask(nx_dst,ny_dst) )
          curr_dst_mask = 1
          curr_dst_mask(1:6,6) = (/ 1, 1, 1, 1, 1, 1 /)
          curr_dst_mask(1:6,5) = (/ 0, 0, 0, 0, 1, 1 /)
          curr_dst_mask(1:6,4) = (/ 1, 0, 0, 1, 1, 1 /)
          curr_dst_mask(1:6,3) = (/ 1, 1, 1, 1, 1, 1 /)
          curr_dst_mask(1:6,2) = (/ 1, 1, 1, 1, 0, 0 /)
          curr_dst_mask(1:6,1) = (/ 1, 1, 1, 1, 0, 0 /)
          call add_grid(curr_dst_grid_id, dst_grid, imask=curr_dst_mask, cartesian=.true.)

          description(1) = "method = conservative"
          description(2) = "Norm type = fracarea"
          strng = " "
          write(strng,'(a,l1)')"Renormalize destarea = ",renorm_destarea
          description(3) = trim(strng)

        case (6)
          member = member + 1
          regrid_method = "bilinear"
          this_norm_opt = "destarea"
          renorm_destarea = .true.
          use_src_mask = .true.
          use_dst_mask = .true.

          !--- Reset the dst mask to its original value
          curr_dst_mask = dst_mask
          call add_grid(curr_dst_grid_id, dst_grid, imask=curr_dst_mask, cartesian=.true.)

          description(1) = "method = bilinear"
          description(2) = "Norm type = destarea"
          strng = " "
          write(strng,'(a,l1)')"Renormalize destarea = ",renorm_destarea
          description(3) = trim(strng)

        case (7)
          member = member + 1
          regrid_method = "patch"
          this_norm_opt = "destarea"
          renorm_destarea = .true.
          use_src_mask = .true.
          use_dst_mask = .true.

          description(1) = "method = patch"
          description(2) = "Norm type = destarea"
          strng = " "
          write(strng,'(a,l1)')"Renormalize destarea = ",renorm_destarea
          description(3) = trim(strng)

        case (8)
          member = member + 1
          regrid_method = "nearest_stod"
          this_norm_opt = "destarea"
          renorm_destarea = .true.
          use_src_mask = .true.
          use_dst_mask = .true.

          description(1) = "method = nearest_stod"
          description(2) = "Norm type = destarea"
          strng = " "
          write(strng,'(a,l1)')"Renormalize destarea = ",renorm_destarea
          description(3) = trim(strng)

        case (9)
          member = member + 1
          regrid_method = "nearest_dtos"
          this_norm_opt = "destarea"
          renorm_destarea = .true.
          use_src_mask = .true.
          use_dst_mask = .true.

          description(1) = "method = nearest_dtos"
          description(2) = "Norm type = destarea"
          strng = " "
          write(strng,'(a,l1)')"Renormalize destarea = ",renorm_destarea
          description(3) = trim(strng)

        case default
          cycle
      end select

      write(6,*)"regrid: member=",member,"  regrid_method=",trim(regrid_method), &
                "  norm_type=",trim(this_norm_opt),"  renormalize_destarea=",renorm_destarea

      !--- Interpolate
      !--- Weights will be calculated if this regrid operation has not yet been seen
      call regrid_Field(src_vinfo%name, src_data, dst_vinfo%name, dst_data, regridID=regrid_id, &
                        src_grid_id=curr_src_grid_id, dst_grid_id=curr_dst_grid_id, &
                        method=regrid_method, norm_type=this_norm_opt, line_type="cartesian", &
                        use_src_mask=use_src_mask, use_dst_mask=use_dst_mask, &
                        src_mask_values=src_mask_values, dst_mask_values=dst_mask_values, &
                        missing_value=0.0_8, renormalize_destarea=renorm_destarea )

      !--- Reformat destination data as a 2D array
      dst_data2d = reshape(dst_data, (/nx_dst,ny_dst/))

      !--- Define an array containing the interpolation error
      !--- determined by comparison with the known function
      if ( associated(err1_data) ) deallocate(err1_data)
      allocate( err1_data(nx_dst,ny_dst) )
      do ix=1,nx_dst
        do iy=1,ny_dst
          err1_data(ix,iy) = test_function( k_test_func, dst_lon(ix), dst_lat(iy) ) - dst_data2d(ix,iy)
        enddo
      enddo

      !--- All masked cells will have zero error
      err1_data = err1_data * curr_dst_mask

      !--- Re-interpolate back onto the source grid and put the result into err2_data
      call regrid_Field(dst_vinfo%name, dst_data, src_vinfo%name, src_data, regridID=err2_regrid_id, &
                        src_grid_id=curr_dst_grid_id, dst_grid_id=curr_src_grid_id, &
                        method=regrid_method, norm_type=this_norm_opt, line_type="cartesian", &
                        use_src_mask=use_dst_mask, use_dst_mask=use_src_mask, &
                        src_mask_values=dst_mask_values, dst_mask_values=src_mask_values, &
                        missing_value=0.0_8, renormalize_destarea=renorm_destarea )

      !--- Reformat the resulting data as a 2D array
      if ( associated(err2_data) ) deallocate(err2_data)
      allocate( err2_data(nx_src,ny_src) )
      err2_data = reshape(src_data, (/nx_src,ny_src/))

      !--- All masked cells will have zero error
      err2_data = err2_data * curr_src_mask

      !--- Put the re-interpolation error into err2_data, defined as the difference
      !--- between the re-interpolated data and the original source data
write(6,*)"regrid:                src_data = ",src_data2d
write(6,*)"regrid: reinterpolated src_data = ",err2_data
      err2_data = src_data2d - err2_data

      !--- Get src_frac and dst_frac from the regrid list create by the call to regrid_field
      if ( associated(regrid_list(regrid_id)%src_frac) ) then
        nx_src_frac = size(regrid_list(regrid_id)%src_frac, dim=1)
        ny_src_frac = size(regrid_list(regrid_id)%src_frac, dim=2)
        if ( associated(src_frac) ) deallocate(src_frac)
        allocate( src_frac(nx_src_frac, ny_src_frac) )
        src_frac = regrid_list(regrid_id)%src_frac
        if ( nx_src_frac == nx_src .and. ny_src_frac == ny_src ) then
          use_src_frac = .true.
        else
          use_src_frac = .false.
          write(6,*)"shape of src_frac differs from shape of src_data."
          write(6,*)"    shape(src_frac) = ",shape(src_frac)
          write(6,*)"    shape(src_data) = ",shape(src_data2d)
        endif
      else
        write(6,*)"regrid: src_frac is not available."
        use_src_frac = .false.
      endif

      if ( associated(regrid_list(regrid_id)%dst_frac) ) then
        nx_dst_frac = size(regrid_list(regrid_id)%dst_frac, dim=1)
        ny_dst_frac = size(regrid_list(regrid_id)%dst_frac, dim=2)
        if ( associated(dst_frac) ) deallocate(dst_frac)
        allocate( dst_frac(nx_dst_frac, ny_dst_frac) )
        dst_frac = regrid_list(regrid_id)%dst_frac
        if ( nx_dst_frac == nx_dst .and. ny_dst_frac == ny_dst ) then
          use_dst_frac = .true.
        else
          use_dst_frac = .false.
          write(6,*)"shape of dst_frac differs from shape of dst_data."
          write(6,*)"    shape(dst_frac) = ",shape(dst_frac)
          write(6,*)"    shape(dst_data) = ",shape(dst_data2d)
        endif
      else
        write(6,*)"regrid: dst_frac is not available."
        use_dst_frac = .false.
      endif

      if ( trim(adjustl(regrid_method)) .ne. "conservative" ) then
        !--- src_frac and dst_frac are not defined
        use_src_frac = .false.
        use_dst_frac = .false.
      endif

      !--- This is the average over non-masked grid cells
      if ( use_src_frac ) then
        src_avg = sum(src_area*src_frac*src_data2d*curr_src_mask)
      else
        src_avg = sum(src_area*src_data2d*curr_src_mask)
      endif
      if ( use_dst_frac ) then
        dst_avg = sum(dst_area*dst_frac*dst_data2d*curr_dst_mask)
      else
        dst_avg = sum(dst_area*dst_data2d*curr_dst_mask)
      endif
      write(6,*)member," src_avg = ",src_avg
      write(6,*)member," dst_avg = ",dst_avg

      !--- Assign certain description lines
      strng = " "
      if ( curr_src_grid_id /= id_src_grid ) then
        !--- This grid mask has been modified from the original
        write(strng,'(a,i2,3a)')"src Grid ID = ",curr_src_grid_id, &
                                "  name = ",trim(regrid_list(regrid_id)%src_grid_name), &
                                "  ..new mask"
      else
        write(strng,'(a,i2,2a)')"src Grid ID = ",curr_src_grid_id, &
                                "  name = ",trim(regrid_list(regrid_id)%src_grid_name)
      endif
      description(4) = trim(strng)

      if ( curr_dst_grid_id /= id_dst_grid ) then
        !--- This grid mask has been modified from the original
        write(strng,'(a,i2,3a)')"dst Grid ID = ",curr_dst_grid_id, &
                                "  name = ",trim(regrid_list(regrid_id)%dst_grid_name), &
                                "  ..new mask"
      else
        write(strng,'(a,i2,2a)')"dst Grid ID = ",curr_dst_grid_id, &
                                "  name = ",trim(regrid_list(regrid_id)%dst_grid_name)
      endif
      description(5) = trim(strng)

      write(strng,'(a,f16.10)')"src Global avg = ",src_avg
      description(6) = trim(strng)
      write(strng,'(a,f16.10)')"dst Global avg = ",dst_avg
      description(7) = trim(strng)

      !--- Add a visual seperator for the last line of each description
      strng = " "
      write(strng,'(a,i2)')"------ END member ",member
      description(8) = trim(strng)

      if ( verbose > 1 ) then
        write(6,*)"src_data: "
        do iy=ny_src,1,-1
          write(6,'(20g18.8)')src_data2d(1:nx_src,iy)
        enddo

        write(6,*)"dst_data: "
        do iy=ny_dst,1,-1
          write(6,'(20g18.8)')dst_data2d(1:nx_dst,iy)
        enddo
      endif

      !--- Write the current description
      call ncdf_write_var(fid, description_id, vals=description, nt=member)

      !--- Write fields
      call ncdf_write_var(fid, src_data_id, vals=src_data2d, nt=member)
      call ncdf_write_var(fid, dst_data_id, vals=dst_data2d, nt=member)

      call ncdf_write_var(fid, err1_data_id, vals=err1_data, nt=member)
      call ncdf_write_var(fid, err2_data_id, vals=err2_data, nt=member)

      call ncdf_write_var(fid, src_mask_id, vals=curr_src_mask, nt=member)
      call ncdf_write_var(fid, dst_mask_id, vals=curr_dst_mask, nt=member)

      if ( use_src_frac ) then
        call ncdf_write_var(fid, src_frac_id, vals=src_frac, nt=member)
      endif
      if ( use_dst_frac ) then
        call ncdf_write_var(fid, dst_frac_id, vals=dst_frac, nt=member)
      endif

    enddo

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

  end subroutine test_grid

  function test_function(kfunc, x, y) result(val)
    integer :: kfunc
    real(8) :: x, y
    real(8) :: val

    !--- Local
    real(8) :: xmin, xmax, ymin, ymax
    real(8) :: xfact, yfact, xn, yn
    real(8) :: x1, x2, y1, y2
    real(8), allocatable :: xb(:), yb(:), fld(:,:)
    real(8) :: pi = 3.141592653589793238_8
    integer :: verbose = 0
    integer :: indx, ix, iy

    xmin = test_func_xmin
    xmax = test_func_xmax
    ymin = test_func_ymin
    ymax = test_func_ymax

    !--- xfact is defined such that xfact*(x-xmin)*(xmax-x) will
    !--- have a maximum value of 1 in the range xmin..xmax
    xfact = 4.0_8/(xmax-xmin)**2

    !--- yfact is defined such that yfact*(y-ymin)*(ymax-y) will
    !--- have a maximum value of 1 in the range ymin..ymax
    yfact = 4.0_8/(ymax-ymin)**2

    select case (kfunc)
      case (1)
        !--- Return a constant value of one
        val = 1.0_8

      case (2)
        !--- Return a 2d quadratic function, constant on all boundaries
        val = 1.0_8 + 0.5_8 * ( ( xfact*(x-xmin)*(xmax-x) ) + &
                                ( yfact*(y-ymin)*(ymax-y) ) )

      case (3)
        !--- Return a 2d trig function (cos*cos) bimodal in x and y
        xn = 1.0_8
        yn = 1.0_8
        val = 1.0_8 + cos( xfact*(x-xmin)*(xmax-x)*xn*pi ) * &
                      cos( yfact*(y-ymin)*(ymax-y)*yn*pi )

      case (4)
        !--- Return a 2d trig function (cos*cos) higher frequency
        xn = 5.0_8
        yn = 5.0_8
        val = 1.0_8 + cos( xfact*(x-xmin)*(xmax-x)*xn*pi ) * &
                      cos( yfact*(y-ymin)*(ymax-y)*yn*pi )

      case (5)
        !--- Return a specific set of values on a 3x3 grid
        !--- Define the x and y cell boundaries
        if ( allocated(xb) ) deallocate(xb)
        allocate( xb(4) )
        if ( allocated(yb) ) deallocate(yb)
        allocate( yb(4) )
        xb(1:4) = (/ 0.0_8, 2.0_8, 4.0_8, 6.0_8 /)
        yb(1:4) = (/ 0.0_8, 2.0_8, 4.0_8, 6.0_8 /)

        !--- Define the x and y cell values
        if ( allocated(fld) ) deallocate(fld)
        allocate( fld(3,3) )
        fld(:,3) = (/  1.0_8, 10.0_8,  1.0_8 /)
        fld(:,2) = (/ 10.0_8,  1.0_8, 10.0_8 /)
        fld(:,1) = (/  1.0_8, 10.0_8,  1.0_8 /)
        ix = 0
        iy = 0
        do idx1=1,size(xb)-1
          if ( x >= xb(idx1) .and. x <=  xb(idx1+1) ) ix = idx1
        enddo
        do idx1=1,size(yb)-1
          if ( y >= yb(idx1) .and. y <=  yb(idx1+1) ) iy = idx1
        enddo
        if ( ix < 1 .or. ix > size(fld,dim=1) .or. iy < 1 .or. iy > size(fld,dim=2) ) then
          write(6,*)"test_function: kfunc=5 ix or iy are out of bounds."
          write(6,*)"  ix=",ix,"  iy=",iy
          call xit("test_function",-1)
        endif

        val = fld(ix,iy)

        deallocate( fld, xb, yb )

        if ( verbose > -1 ) then
          write(6,*)"test_function: x=",x,"  y=",y,"  ix=",ix,"  iy=",iy,"  val=",val
        endif

      case (6)
        !--- Return a specific set of values on a 5x5 grid
        !--- Define the x and y cell boundaries
        if ( allocated(xb) ) deallocate(xb)
        allocate(xb(6))
        if ( allocated(yb) ) deallocate(yb)
        allocate(yb(6))
        xb(1:6) = (/ 0.0_8, 1.2_8, 2.4_8, 3.6_8, 4.8_8, 6.0_8 /)
        yb(1:6) = (/ 0.0_8, 1.2_8, 2.4_8, 3.6_8, 4.8_8, 6.0_8 /)

        !--- Define the x and y cell values
        if ( allocated(fld) ) deallocate(fld)
        allocate( fld(5,5) )
        fld(:,5) = (/  1.0_8, 10.0_8,  1.0_8, 10.0_8,  1.0_8 /)
        fld(:,4) = (/ 10.0_8,  1.0_8, 10.0_8,  1.0_8, 10.0_8 /)
        fld(:,3) = (/  1.0_8, 10.0_8,  1.0_8, 10.0_8,  1.0_8 /)
        fld(:,2) = (/ 10.0_8,  1.0_8, 10.0_8,  1.0_8, 10.0_8 /)
        fld(:,1) = (/  1.0_8, 10.0_8,  1.0_8, 10.0_8,  1.0_8 /)
        ix = 0
        iy = 0
        do idx1=1,size(xb)-1
          if ( x >= xb(idx1) .and. x <=  xb(idx1+1) ) ix = idx1
        enddo
        do idx1=1,size(yb)-1
          if ( y >= yb(idx1) .and. y <=  yb(idx1+1) ) iy = idx1
        enddo
        if ( ix < 1 .or. ix > size(fld,dim=1) .or. iy < 1 .or. iy > size(fld,dim=2) ) then
          write(6,*)"test_function: kfunc=5 ix or iy are out of bounds."
          write(6,*)"  ix=",ix,"  iy=",iy
          call xit("test_function",-1)
        endif

        val = fld(ix,iy)

        deallocate( fld, xb, yb )

        if ( verbose > -1 ) then
          write(6,*)"test_function: x=",x,"  y=",y,"  ix=",ix,"  iy=",iy,"  val=",val
        endif

      case (7)
        !--- Return a specific set of values on a 6x6 grid
        !--- Define the x and y cell boundaries
        if ( allocated(xb) ) deallocate(xb)
        allocate(xb(7))
        if ( allocated(yb) ) deallocate(yb)
        allocate(yb(7))
        xb(1:7) = (/ 0.0_8, 1.0_8, 2.0_8, 3.0_8, 4.0_8, 5.0_8, 6.0_8 /)
        yb(1:7) = (/ 0.0_8, 1.0_8, 2.0_8, 3.0_8, 4.0_8, 5.0_8, 6.0_8 /)

        !--- Define the x and y cell values
        if ( allocated(fld) ) deallocate(fld)
        allocate( fld(6,6) )

        !--- This is an example of a land mask
        fld(:,6) = (/ 1.0_8, 1.0_8, 1.0_8, 1.0_8, 1.0_8, 1.0_8 /)
        fld(:,5) = (/ 0.0_8, 0.0_8, 0.0_8, 0.0_8, 1.0_8, 1.0_8 /)
        fld(:,4) = (/ 1.0_8, 0.0_8, 0.0_8, 1.0_8, 1.0_8, 1.0_8 /)
        fld(:,3) = (/ 1.0_8, 1.0_8, 1.0_8, 1.0_8, 1.0_8, 1.0_8 /)
        fld(:,2) = (/ 1.0_8, 1.0_8, 1.0_8, 1.0_8, 0.0_8, 0.0_8 /)
        fld(:,1) = (/ 1.0_8, 1.0_8, 1.0_8, 1.0_8, 0.0_8, 0.0_8 /)

        ix = 0
        iy = 0
        do idx1=1,size(xb)-1
          if ( x >= xb(idx1) .and. x <=  xb(idx1+1) ) ix = idx1
        enddo
        do idx1=1,size(yb)-1
          if ( y >= yb(idx1) .and. y <=  yb(idx1+1) ) iy = idx1
        enddo
        if ( ix < 1 .or. ix > size(fld,dim=1) .or. iy < 1 .or. iy > size(fld,dim=2) ) then
          write(6,*)"test_function: kfunc=7 ix or iy are out of bounds."
          write(6,*)"  ix=",ix,"  iy=",iy
          call xit("test_function",-1)
        endif

        val = fld(ix,iy)

        deallocate( fld, xb, yb )

        if ( verbose > -1 ) then
          write(6,*)"test_function: x=",x,"  y=",y,"  ix=",ix,"  iy=",iy,"  val=",val
        endif

      case (8)
        !--- Return a specific set of values on a 6x6 grid
        !--- Define the x and y cell boundaries
        if ( allocated(xb) ) deallocate(xb)
        allocate(xb(7))
        if ( allocated(yb) ) deallocate(yb)
        allocate(yb(7))
        xb(1:7) = (/ 0.0_8, 1.0_8, 2.0_8, 3.0_8, 4.0_8, 5.0_8, 6.0_8 /)
        yb(1:7) = (/ 0.0_8, 1.0_8, 2.0_8, 3.0_8, 4.0_8, 5.0_8, 6.0_8 /)

        !--- Define the x and y cell values
        if ( allocated(fld) ) deallocate(fld)
        allocate( fld(6,6) )

        !--- This is an example of SICN on a 6x6 grid
        fld(:,6) = (/ 0.5_8, 0.5_8, 1.0_8,  0.0_8, 0.0_8, 1.0_8 /)
        fld(:,5) = (/ 0.0_8, 0.0_8, 0.0_8,  0.0_8, 1.0_8, 0.0_8 /)
        fld(:,4) = (/ 0.0_8, 0.0_8, 0.0_8, 0.15_8, 0.1_8, 0.1_8 /)
        fld(:,3) = (/ 0.3_8, 0.0_8, 0.0_8, 0.15_8, 0.1_8, 0.1_8 /)
        fld(:,2) = (/ 0.0_8, 0.0_8, 0.0_8,  0.3_8, 0.0_8, 0.0_8 /)
        fld(:,1) = (/ 0.0_8, 0.0_8, 0.0_8,  0.1_8, 0.0_8, 0.0_8 /)

        ix = 0
        iy = 0
        do idx1=1,size(xb)-1
          if ( x >= xb(idx1) .and. x <=  xb(idx1+1) ) ix = idx1
        enddo
        do idx1=1,size(yb)-1
          if ( y >= yb(idx1) .and. y <=  yb(idx1+1) ) iy = idx1
        enddo
        if ( ix < 1 .or. ix > size(fld,dim=1) .or. iy < 1 .or. iy > size(fld,dim=2) ) then
          write(6,*)"test_function: kfunc=8 ix or iy are out of bounds."
          write(6,*)"  ix=",ix,"  iy=",iy
          call xit("test_function",-1)
        endif

        val = fld(ix,iy)

        deallocate( fld, xb, yb )

        if ( verbose > -1 ) then
          write(6,*)"test_function: x=",x,"  y=",y,"  ix=",ix,"  iy=",iy,"  val=",val
        endif

      case default
        write(6,*)"test_function: Invalid function kind = ",kfunc
        call xit("test_function",-9)

    end select

    if ( verbose > 0 ) then
      write(6,*)"test_function: x=",x,"  y=",y,"  val=",val
      if ( verbose > 1 ) then
        write(6,*)"test_function: xfact=",xfact,"  yfact=",yfact
        write(6,*)"test_function: argx/(n*pi)=",xfact*(x-xmin)*(xmax-x), &
                               "  argy/(n*pi)=",yfact*(y-ymin)*(ymax-y)
        write(6,*)"test_function:  sin(argx)=",sin(xfact*(x-xmin)*(xmax-x)*pi), &
                                "  sin(argy)=",sin(yfact*(y-ymin)*(ymax-y)*pi)
      endif
    endif

  end function test_function

end module regrid_mod

program regrid

  !--- Interpolate data read from input files onto another grid
  !--- Input and output grids are specified by the user

  !--- ccc_esmf also uses com_cpl, strings, ccc_netcdf, ccc_io, store
  use regrid_mod

  implicit none

  character(64) :: grid_id
  logical :: ok, src_overlap, dst_overlap
  integer :: src_dims(2)

  integer :: nc4to8
  integer :: ios
  integer, parameter :: iuin=1, iuout=2, iuorig=3
  logical :: remote, remote_src_grid, remote_dst_grid
  logical :: match_found, testing=.false.
  character(4) :: cckind, ccname

  real(kind=8), allocatable :: fld_in(:,:), fld_out(:,:)
  real(kind=8), allocatable :: src_arr(:), dst_arr(:), src_arr2(:), dst_arr2(:)

  !--- Index of the current grid description returned from create_grid_desc
  integer :: curr_grid

  character(256) :: area_frac_file

  type(ESMF_Field) :: src_area_field, dst_area_field

  !--- CCCma io specific
  integer :: machine, me32o64
  integer :: intsize, inteflt
  external me32o64, inteflt
  common /machtyp/ machine,intsize

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

  if (use_timers) call timer_cpu( 'start', idx=1, tag='regrid' )

  ! Initialize machine and intsize (required for CCCma io routines)
  machine = me32o64(nx)
  intsize = inteflt(nx)

  !--- Output file name
  fout = 'regrid_output'
  usr_def_fout = .false.

  !--- File name to hold original data selected from input data
  forig = 'regrid_orig'

  !--- Files containing descriptions of the source and destination grids
  src_grid=" "
  dst_grid=" "

  !--- The name of a CCCma binary file containing ocean data description info
  ocn_desc = " "

  !--- The name of a CCCma binary file containing a gaussian grid land mask
  !---   -1 for land, anything else is not-land
  mask_gg = " "

  !--- The name of a netcdf file (used by NEMO) containing ORCA coordinates
  nemo_coords = " "

  !--- (start|stop)_ibuf2, if >0, will indicate the first and last value of ibuf2
  !--- to be read from the input file. That is, any records with ibuf2 in the range
  !--- start_ibuf2 <= ibuf2 <= stop_ibuf2 will be kept, all otehr will be ignored
  start_ibuf2 = -1
  stop_ibuf2  = -1

  !--- Process the command line args
  call process_clg()

  if ( create_frac_land_mask ) then
    !--- Set defaults for src/dst grids when creating fractional land mask
    if ( len_trim(src_grid) < 1 ) then
      !--- Define src_grid here if not supplied on the command line
      src_grid = "t63"
    endif
    if ( len_trim(dst_grid) < 1 ) then
      !--- Define dst_grid here if not supplied on the command line
      dst_grid = "orca1"
    endif
  endif

  testing = .false.
  if ( src_grid(1:4) .eq. "test" ) testing = .true.
  if ( dst_grid(1:4) .eq. "test" ) testing = .true.

  select case ( trim(adjustl(src_grid)) )
    case ( 'atm' )
      !--- Use a standard atm grid desciption file
      src_grid = "grid_desc_t63_o_gcm3.5_landmask_128x64v5.nc"
    case ( 'ocn', 'nemo' )
      !--- Use a standard ocean grid description file
      src_grid = "grid_desc_orca1_nemo_3.4_orca1_coordinates_ukorca1_with_mask.nc"
  end select

  select case ( trim(adjustl(dst_grid)) )
    case ( 'atm' )
      !--- Use a standard atm grid desciption file
      dst_grid = "grid_desc_t63_o_gcm3.5_landmask_128x64v5.nc"
    case ( 'ocn', 'nemo' )
      !--- Use a standard ocean grid description file
      dst_grid = "grid_desc_orca1_nemo_3.4_orca1_coordinates_ukorca1_with_mask.nc"
  end select

  if ( grid_desc_only ) then
    !--- Create a grid description file then exit
    if ( len_trim(src_grid) .le. 0 ) then
      write(6,*)'src_grid must be defined when only grid description file is created.'
      call usage(-4)
    endif

  else
    !--- Normal processing
    if ( .not. create_frac_land_mask .and. .not. testing .and. nfiles .le. 0) then
      write(6,*)'regrid: No input data files were specified on the command line.'
      call usage(-4)
    endif

    if ( len_trim(src_grid) .le. 0 ) then
      ok = .false.
      if ( len_trim(dst_grid) > 0 .and. len_trim(ocn_desc) > 0 ) then
        !--- If the destination grid is provided along with an ocean description file
        !--- then assume the source grid is described by ocn_desc
        src_grid = 'use_ocn_desc'
        ok = .true.
      endif
      if ( .not. ok ) then
        write(6,*)'regrid: No source grid file was specified on the command line.'
        call usage(-4)
      endif
    endif

    if ( len_trim(dst_grid) .le. 0 ) then
      ok = .false.
      if ( len_trim(src_grid) > 0 .and. len_trim(ocn_desc) > 0 ) then
        !--- If the source grid is provided along with an ocean description file
        !--- then assume the destination grid is described by ocn_desc
        dst_grid = 'use_ocn_desc'
        ok = .true.
      endif
      if ( .not. ok ) then
        write(6,*)'regrid: No destination grid file was specified on the command line.'
        call usage(-4)
      endif
    endif
  endif

  !--- Verify that we have a valid map type
  select case(trim(regrid_method))
    case ('conservative', 'cv')
      regrid_method = 'conservative'
      if (.not.usr_def_fout) fout="regrid_output_cv"
    case ('bilinear', 'bl')
      regrid_method = 'bilinear'
      if (.not.usr_def_fout) fout="regrid_output_bl"
    case ('patch')
      regrid_method = 'patch'
      if (.not.usr_def_fout) fout="regrid_output_patch"
    case ('nearest_stod', 'stod')
      regrid_method = 'nearest_stod'
      if (.not.usr_def_fout) fout="regrid_output_stod"
    case ('nearest_dtos', 'dtos')
      regrid_method = 'nearest_dtos'
      if (.not.usr_def_fout) fout="regrid_output_dtos"
    case default
      write(6,'(2a)')'regrid: Invalid map_type ',trim(regrid_method)
      call usage(-4)
  end select

  !--- Verify that we have a valid normalization option
  select case(trim(this_norm_opt))
    case ('fracarea')
    case ('destarea')
    case ('none')
    case default
      write(6,'(2a)')'regrid: Invalid norm_opt ',trim(this_norm_opt)
      call usage(-4)
  end select

  !--- Valid grid IDs that are recognized by create_grid_desc are:
  !---    t31, t47, t63, t95, t127
  !---    ocnt_256x192, ocnv_256x192, use_ocn_desc
  !---    orca2, orca1, orca025

  !--- Ensure the input grid description files exist, if not then create it
  if ( .not. locate_file(trim(src_grid),remote_src_grid) ) then
    grid_id = trim(adjustl(src_grid))
    if ( .not. src_masked ) then
      curr_mask_type = -1
    else
      curr_mask_type = 0
    endif
    !--- Create a grid description file based on the value of grid_id
    !--- create_grid_desc will return a new name for src_grid
    if ( len_trim(ocn_desc) > 0 .and. &
         len_trim(mask_gg) > 0 .and. &
         len_trim(nemo_coords) > 0 ) then
      !--- Pass nemo_coords, ocn_desc and mask_gg even though only one will be used
      !--- depending on the value of grid_id
      call create_grid_desc( grid_id, src_grid, curr_grid, mask_type=curr_mask_type, &
            ocn_desc=trim(ocn_desc), mask_gg=trim(mask_gg), nemo_coords=trim(nemo_coords), &
            olaplon=olaplon )
    else if ( len_trim(ocn_desc) > 0 .and. len_trim(mask_gg) > 0 ) then
      !--- Pass both ocn_desc and mask_gg even though only one will be used
      !--- depending on the value of grid_id
      call create_grid_desc( grid_id, src_grid, curr_grid, mask_type=curr_mask_type, &
                             ocn_desc=trim(ocn_desc), mask_gg=trim(mask_gg), &
                             olaplon=olaplon )
    else if ( len_trim(ocn_desc) > 0 .and. len_trim(nemo_coords) > 0 ) then
      !--- Pass both ocn_desc and nemo_coords even though only one will be used
      !--- depending on the value of grid_id
      call create_grid_desc( grid_id, src_grid, curr_grid, mask_type=curr_mask_type, &
                             ocn_desc=trim(ocn_desc), nemo_coords=trim(nemo_coords), &
                             olaplon=olaplon )
    else if ( len_trim(mask_gg) > 0 .and. len_trim(nemo_coords) > 0 ) then
      !--- Pass both mask_gg and nemo_coords even though only one will be used
      !--- depending on the value of grid_id
      call create_grid_desc( grid_id, src_grid, curr_grid, mask_type=curr_mask_type, &
                             mask_gg=trim(mask_gg), nemo_coords=trim(nemo_coords), &
                             olaplon=olaplon )
    else if ( len_trim(mask_gg) > 0 ) then
      call create_grid_desc( grid_id, src_grid, curr_grid, mask_type=curr_mask_type, &
                             mask_gg=trim(mask_gg), olaplon=olaplon )
    else if ( len_trim(ocn_desc) > 0 ) then
      call create_grid_desc( grid_id, src_grid, curr_grid, mask_type=curr_mask_type, &
                             ocn_desc=trim(ocn_desc), olaplon=olaplon )
    else if ( len_trim(nemo_coords) > 0 ) then
      call create_grid_desc( grid_id, src_grid, curr_grid, mask_type=curr_mask_type, &
                             nemo_coords=trim(nemo_coords), olaplon=olaplon )
    else
      call create_grid_desc( grid_id, src_grid, curr_grid, mask_type=curr_mask_type, &
                             olaplon=olaplon )
    endif
  endif

  if ( grid_desc_only ) then
    if ( remote_src_grid ) then
      strng=' '
      strng="release "//trim(src_grid)
      call system(trim(strng))
    endif
    stop
  endif

  !--- Ensure the output grid description files exist, if not then create it
  if ( .not. locate_file(trim(dst_grid),remote_dst_grid) ) then
    grid_id = trim(adjustl(dst_grid))
    if ( .not. dst_masked ) then
      curr_mask_type = -1
    else
      curr_mask_type = 0
    endif
    !--- Create a grid description file based on the value of grid_id
    !--- create_grid_desc will return a new name for dst_grid
    if ( len_trim(ocn_desc) > 0 .and. &
         len_trim(mask_gg) > 0 .and. &
         len_trim(nemo_coords) > 0 ) then
      !--- Pass nemo_coords, ocn_desc and mask_gg even though only one will be used
      !--- depending on the value of grid_id
      call create_grid_desc( grid_id, dst_grid, curr_grid, mask_type=curr_mask_type, &
            ocn_desc=trim(ocn_desc), mask_gg=trim(mask_gg), nemo_coords=trim(nemo_coords), &
            olaplon=olaplon )
    else if ( len_trim(ocn_desc) > 0 .and. len_trim(mask_gg) > 0 ) then
      !--- Pass both ocn_desc and mask_gg even though only one will be used
      !--- depending on the value of grid_id
      call create_grid_desc( grid_id, dst_grid, curr_grid, mask_type=curr_mask_type, &
                             ocn_desc=trim(ocn_desc), mask_gg=trim(mask_gg), &
                             olaplon=olaplon )
    else if ( len_trim(ocn_desc) > 0 .and. len_trim(nemo_coords) > 0 ) then
      !--- Pass both ocn_desc and nemo_coords even though only one will be used
      !--- depending on the value of grid_id
      call create_grid_desc( grid_id, dst_grid, curr_grid, mask_type=curr_mask_type, &
                             ocn_desc=trim(ocn_desc), nemo_coords=trim(nemo_coords), &
                             olaplon=olaplon )
    else if ( len_trim(mask_gg) > 0 .and. len_trim(nemo_coords) > 0 ) then
      !--- Pass both mask_gg and nemo_coords even though only one will be used
      !--- depending on the value of grid_id
      call create_grid_desc( grid_id, dst_grid, curr_grid, mask_type=curr_mask_type, &
                             mask_gg=trim(mask_gg), nemo_coords=trim(nemo_coords), &
                             olaplon=olaplon )
    else if ( len_trim(mask_gg) > 0 ) then
      call create_grid_desc( grid_id, dst_grid, curr_grid, mask_type=curr_mask_type, &
                             mask_gg=trim(mask_gg), olaplon=olaplon )
    else if ( len_trim(ocn_desc) > 0 ) then
      call create_grid_desc( grid_id, dst_grid, curr_grid, mask_type=curr_mask_type, &
                             ocn_desc=trim(ocn_desc), olaplon=olaplon )
    else if ( len_trim(nemo_coords) > 0 ) then
      call create_grid_desc( grid_id, dst_grid, curr_grid, mask_type=curr_mask_type, &
                             nemo_coords=trim(nemo_coords), olaplon=olaplon )
    else
      call create_grid_desc( grid_id, dst_grid, curr_grid, mask_type=curr_mask_type, &
                               olaplon=olaplon )
    endif
  endif

  !--- Initialize ESMF
  call ESMF_Initialize ( defaultCalKind=ESMF_CALKIND_GREGORIAN, &
                          defaultlogfilename="Regrid.Log", rc=rc)
  if (rc /= ESMF_SUCCESS) then
    write(6,*)'regrid: Problem initializing ESMF.'
    call ESMF_Finalize(endflag=ESMF_END_ABORT)
    call xit("regrid",-1)
  endif

  if ( create_frac_land_mask ) then
    call create_fractional_land_mask()

    !--- Remove any links to the grid definition files
    if ( remote_src_grid ) then
      strng=' '
      strng="release "//trim(src_grid)
      call system(trim(strng))
    endif
    if ( remote_dst_grid ) then
      strng=' '
      strng="release "//trim(dst_grid)
      call system(trim(strng))
    endif

    stop
  endif

  if ( testing ) then
    write(6,*)"About to test things"
    call test_grid()
    stop
  endif

  !--- Add a definition for the source grid (internal to ccc_esmf)
  call add_grid(id_src_grid, src_grid, &
                grid_dims=src_grid_dims, &
                grid_edges=src_grid_edges, &
                grid_size=src_grid_size)

  !--- Retrieve a pointer to the src grid mask
  !--- Note that reassigning src_grid_mask here will change the mask that is used in
  !--- the remap calculation below, since src_grid_mask is a pointer to that mask
  call ESMF_GridGetItem(grid_list(id_src_grid), &
                        staggerloc=ESMF_STAGGERLOC_CENTER, &
                        itemflag=ESMF_GRIDITEM_MASK, &
                        farrayPtr=src_grid_mask, rc=rc)
  if (rc /= ESMF_SUCCESS) then
    write(6,*)'regrid: Unable to retrive a pointer to the source grid mask.'
    call ESMF_Finalize(endflag=ESMF_END_ABORT)
    call xit("regrid",-2)
  endif
  if ( .not. associated(src_grid_mask) ) then
    write(6,*)'regrid: Pointer to the source grid mask is not associated.'
    call ESMF_Finalize(endflag=ESMF_END_ABORT)
    call xit("regrid",-2)
  endif
  if ( size(src_grid_mask,1) /= src_grid_dims(1) ) then
    write(6,*)"regrid: First dim of src grid does not match src_grid_dims."
    write(6,*)"   mask_dim(1),src_grid_dims(1) = ",size(src_grid_mask,1),src_grid_dims(1)
    call xit("regrid",-1)
  endif
  if ( size(src_grid_mask,2) /= src_grid_dims(2) ) then
    write(6,*)"regrid: Second dim of src grid does not match src_grid_dims."
    write(6,*)"   mask_dim(2),src_grid_dims(2) = ",size(src_grid_mask,2),src_grid_dims(2)
    call xit("regrid",-1)
  endif

!xxx  !--- src_area_field will be an ESMF field type to contain source field values
!xxx  src_area_field = ESMF_FieldCreate( grid_list(id_src_grid), ESMF_TYPEKIND_R8, rc=rc)
!xxx  if (rc /= ESMF_SUCCESS) then
!xxx    write(6,*)'regrid: Unable to create a source area field.'
!xxx    call ESMF_Finalize(endflag=ESMF_END_ABORT)
!xxx    call xit("regrid",-2)
!xxx  endif
!xxx
!xxx  !--- Compute src grid cell areas and store them in the src_area_field object
!xxx  call ESMF_FieldRegridGetArea(src_area_field, rc=rc)
!xxx  if (rc /= ESMF_SUCCESS) then
!xxx    write(6,*)'regrid: Unable to calculate areas in source field.'
!xxx    call ESMF_Finalize(endflag=ESMF_END_ABORT)
!xxx    call xit("regrid",-2)
!xxx  endif
!xxx
!xxx  !--- Extract the 2D area array from src_area_field
!xxx  if ( associated(src_area_data) ) nullify( src_area_data )
!xxx  call ESMF_FieldGet( src_area_field, farrayPtr=src_area_data, rc=rc )
!xxx  if (rc /= ESMF_SUCCESS) then
!xxx    write(6,*)'regrid: Unable to extract a pointer to area data from source field.'
!xxx    call ESMF_Finalize(endflag=ESMF_END_ABORT)
!xxx    call xit("regrid",-2)
!xxx  endif
!xxx
!xxx  write(6,*)"sum(src_area_data) = ",sum(src_area_data)
!xxx
!xxxwrite(6,*)"regrid: src grid_edges = ",src_grid_edges
!xxxcall flush(6)

  !--- Add a definition for the destination grid (internal to ccc_esmf)
  call add_grid(id_dst_grid, dst_grid, &
                grid_dims=dst_grid_dims, &
                grid_edges=dst_grid_edges, &
                grid_size=dst_grid_size)

  !--- Retrieve a pointer to the dst grid mask
  !--- Note that reassigning dst_grid_mask here will change the mask that is used in
  !--- the remap calculation below, since dst_grid_mask is a pointer to that mask
  call ESMF_GridGetItem(grid_list(id_dst_grid), &
                        staggerloc=ESMF_STAGGERLOC_CENTER, &
                        itemflag=ESMF_GRIDITEM_MASK, &
                        farrayPtr=dst_grid_mask, rc=rc)
  if (rc /= ESMF_SUCCESS) then
    write(6,*)'regrid: Unable to retrive a pointer to the source grid mask.'
    call ESMF_Finalize(endflag=ESMF_END_ABORT)
    call xit("regrid",-2)
  endif
  if ( .not. associated(dst_grid_mask) ) then
    write(6,*)'regrid: Pointer to the source grid mask is not associated.'
    call ESMF_Finalize(endflag=ESMF_END_ABORT)
    call xit("regrid",-2)
  endif
  if ( size(dst_grid_mask,1) /= dst_grid_dims(1) ) then
    write(6,*)"regrid: First dim of dst grid does not match dst_grid_dims."
    write(6,*)"   mask_dim(1),dst_grid_dims(1) = ",size(dst_grid_mask,1),dst_grid_dims(1)
    call xit("regrid",-1)
  endif
  if ( size(dst_grid_mask,2) /= dst_grid_dims(2) ) then
    write(6,*)"regrid: Second dim of dst grid does not match dst_grid_dims."
    write(6,*)"   mask_dim(2),dst_grid_dims(2) = ",size(dst_grid_mask,2),dst_grid_dims(2)
    call xit("regrid",-1)
  endif

!xxx  !--- dst_area_field will be an ESMF field type to contain destination cell area
!xxx  dst_area_field = ESMF_FieldCreate( grid_list(id_dst_grid), ESMF_TYPEKIND_R8, rc=rc)
!xxx  if (rc /= ESMF_SUCCESS) then
!xxx    write(6,*)'regrid: Unable to create a destination area field.'
!xxx    call ESMF_Finalize(endflag=ESMF_END_ABORT)
!xxx    call xit("regrid",-2)
!xxx  endif
!xxx
!xxx  !--- Compute dst grid cell areas and store them in the dst_area_field object
!xxx  call ESMF_FieldRegridGetArea(dst_area_field, rc=rc)
!xxx  if (rc /= ESMF_SUCCESS) then
!xxx    write(6,*)'regrid: Unable to calculate areas in source field.'
!xxx    call ESMF_Finalize(endflag=ESMF_END_ABORT)
!xxx    call xit("regrid",-2)
!xxx  endif
!xxx
!xxx  !--- Extract the 2D area array from dst_area_field
!xxx  if ( associated(dst_area_data) ) nullify( dst_area_data )
!xxx  call ESMF_FieldGet( dst_area_field, farrayPtr=dst_area_data, rc=rc )
!xxx  if (rc /= ESMF_SUCCESS) then
!xxx    write(6,*)'regrid: Unable to extract a pointer to area data from destination field.'
!xxx    call ESMF_Finalize(endflag=ESMF_END_ABORT)
!xxx    call xit("regrid",-2)
!xxx  endif
!xxx
!xxx  write(6,*)"sum(dst_area_data) = ",sum(dst_area_data)
!xxx
!xxxwrite(6,*)"regrid: dst grid_edges = ",dst_grid_edges
!xxxcall flush(6)

  !--- Determine which grid masks are to be used
  use_src_mask=.false.
  use_dst_mask=.true.
  if ( associated(src_mask_values) ) deallocate(src_mask_values)
  allocate( src_mask_values(1) )
  src_mask_values = 0
  if ( associated(dst_mask_values) ) deallocate(dst_mask_values)
  allocate( dst_mask_values(1) )
  dst_mask_values = 0

  !--- This invocation of find_regrid_id will cause weights to be generated
  !--- for this remapping operation. These weights will be written to a file
  !--- as well as being stored internal to ccc_esmf
  call find_regrid_id(regrid_id, src_grid, dst_grid, regrid_method,  &
                      use_src_mask, use_dst_mask, &
                      src_mask_values, dst_mask_values, &
                      norm_type=this_norm_opt, create=.true.)

  if ( .false. .and. trim(adjustl(regrid_method)) .eq. "conservative" ) then
    !--- This is a conservative remap
    area_frac_file="cell_area_and_fraction"

    !--- Generating regid weights for a conservative remap will also cause
    !--- src and dst cell area and fraction fields to be stored internally
    strng = " "
    write(strng,'(a,i3.3)')"src_frac",regrid_id
    nwrds = store_ReadRecord(trim(strng), other_gid, src_frac, nx=nx_sfrac, ny=ny_sfrac, idx=-1 )
    if ( verbose > 10 ) call store_ListRecord(trim(strng), other_gid)

    strng = " "
    write(strng,'(a,i3.3)')"dst_frac",regrid_id
    nwrds = store_ReadRecord(trim(strng), other_gid, dst_frac, nx=nx_dfrac, ny=ny_dfrac, idx=-1 )
    if ( verbose > 10 ) call store_ListRecord(trim(strng), other_gid)

    !--- Write src and dst cell area and fraction fields to a CCCma binary file
    ibuf(1:8) = 1
    ibuf(1) = transfer("GRID", 1)
    ibuf(7) = 0
    open(701, file=trim(area_frac_file), form='unformatted')
      ibuf(3) = transfer("SRCF", 1)
      ibuf(5) = nx_sfrac
      ibuf(6) = ny_sfrac
      wrk(1:ibuf(5)*ibuf(6)) = src_frac(1:ibuf(5)*ibuf(6))
      where ( wrk(1:ibuf(5)*ibuf(6)) > 1.0 ) wrk = 1.0

if ( .false. ) then
      if ( src_grid_dims(1) /= nx_sfrac .or. src_grid_dims(2) /= ny_sfrac ) then
        write(6,*)"regrid: Size mismatch between internal src_fracNNN and src grid"
        call xit("regrid",-1)
      endif
      allocate( itmp_data(src_grid_dims(1)*src_grid_dims(2)) )
      itmp_data = reshape( src_grid_mask, (/src_grid_dims(1)*src_grid_dims(2)/) )
      where ( itmp_data == 1 ) wrk(1:ibuf(5)*ibuf(6)) = 0.0
endif

      call putfld2(701,wrk,ibuf,maxx)

    close(701)

  endif

  !--- Define cpl_vinfo_t variables (internal to ccc_esmf)
  allocate(src_vinfo)
  src_vinfo%name     = "src_field"
  src_vinfo%grid     = trim(src_grid)
  src_vinfo%tag      = cpl_var_new_tag(5000)
  src_vinfo%size     = src_grid_dims(1)*src_grid_dims(2)
  src_vinfo%nlon     = src_grid_dims(1)  !--- Do not include overlap
  src_vinfo%nlat     = src_grid_dims(2)
  src_vinfo%olap     = 0
  src_vinfo%ncopy    = 1
  src_vinfo%freq     = 0
  src_vinfo%ccc_name = " SRC"
  src_vinfo%regrid%method   = trim(regrid_method)
  src_vinfo%regrid%norm_opt = trim(this_norm_opt)
  src_vinfo%regrid%masked   = use_src_mask
  !--- Add this variable defintion, overwriting any existing def with the same name
  call add_cpl_var(src_vinfo, exists=1)

  allocate(dst_vinfo)
  dst_vinfo%name     = "dst_field"
  dst_vinfo%grid     = trim(dst_grid)
  dst_vinfo%tag      = cpl_var_new_tag(5000)
  dst_vinfo%size     = dst_grid_dims(1)*dst_grid_dims(2)
  dst_vinfo%nlon     = dst_grid_dims(1)  !--- Do not include overlap
  dst_vinfo%nlat     = dst_grid_dims(2)
  dst_vinfo%olap     = 0
  dst_vinfo%ncopy    = 1
  dst_vinfo%freq     = 0
  dst_vinfo%ccc_name = " DST"
  dst_vinfo%regrid%method   = trim(regrid_method)
  dst_vinfo%regrid%norm_opt = trim(this_norm_opt)
  dst_vinfo%regrid%masked   = use_dst_mask
  !--- Add this variable defintion, overwriting any existing def with the same name
  call add_cpl_var(dst_vinfo, exists=1)

  !--- Assign dimension sizes using the grid_dims returned from the subroutine add_grid
  nx1 = src_grid_dims(1)
  ny1 = src_grid_dims(2)
  nx2 = dst_grid_dims(1)
  ny2 = dst_grid_dims(2)

  !--- Attach a new CCCma binary file to hold output data
  if ( .not. nc_out ) then
    open(iuout, file=trim(fout), form='unformatted')
  endif

  if ( write_orig_fld ) then
    !--- Open a separate file to hold the original data
    open(iuorig, file=trim(forig), form='unformatted')
  endif

  !--- Loop over input files and read/interpolate every GRID record found
  READ_FILES: do kfile=1,nfiles
    curr_file = files(kfile)
    call cc_open(iuin,trim(curr_file),remote,1)

    ok = .true.
    nrec = 0
    READ_RECORDS: do while (ok)
      !--- Read the next GRID record from the input file
      call getfld2(iuin, wrk, nc4to8("GRID"), -1, -1, -1, ibuf, maxx, ok)
      if (verbose > 1) then
        write(6,'(a,a4,1x,i10,1x,a4,1x,5i8)')'read: ',ibuf(1:8)
      endif

      if ( .not. ok ) exit READ_RECORDS

      !--- Get KIND and NAME as char from the integer label
      cckind = transfer(ibuf(1),cckind)
      ccname = transfer(ibuf(3),ccname)

      if (nvlist.gt.0) then
        !--- Include only names found in vlist
        match_found = .false.
        do idx1=1,nvlist
          if ( trim(adjustl(ccname)) .eq. &
               trim(adjustl(vlist(idx1))) ) then
            match_found = .true.
            exit
          endif
        enddo
        !--- If this variable name is not in the list then skip
        if ( .not.match_found ) cycle READ_RECORDS
      endif

      !--- Exclude records outside of any user supplied time range
      if ( start_ibuf2 > -1 ) then
        if ( ibuf(2) < start_ibuf2 ) cycle READ_RECORDS
      endif
      if ( stop_ibuf2 > -1 ) then
        if ( ibuf(2) >  stop_ibuf2 ) cycle READ_RECORDS
      endif

      !--- Determine size of input record and compare with input grid size
      nx = ibuf(5)
      ny = ibuf(6)
      if ( (nx /= nx1 .and. nx /= nx1+1) .or. ny /= ny1 ) then
        !--- The current record is not defined on the input grid
        if (verbose > 0) then
          write(6,'(a,a4,1x,i10,1x,a4,1x,5i8)')'Skipping: ',ibuf(1:8)
          write(6,'(a)')'  which does not conform to the input grid size.'
          write(6,'(a,2i8)')'    Input grid nx, ny = ',nx1,ny1
          write(6,'(a,2i8)')'Current record nx, ny = ',nx,ny
        endif
        cycle READ_RECORDS
      endif

      !--- Count the number of GRID records that are read
      nrec = nrec + 1

      if (verbose > 0) then
        write(6,'(i6,a,a4,1x,i10,1x,a4,1x,5i8)')nrec,' interpolating: ',ibuf(1:8)
      endif

!xxx      !--- Allocate and assign the input and output arrays for interpolation
!xxx       if ( allocated(src_arr) ) deallocate(src_arr)
!xxx       if ( allocated(dst_arr) ) deallocate(dst_arr)
!xxx       allocate( src_arr(nx*ny), dst_arr(nx2*ny2) )
!xxx       src_arr(1:nx*ny) = wrk(1:nx*ny)
!xxx       dst_arr(:) = 0.0
!xxx       src_dims(1:2) = (/nx,ny/)

      !--- Interpolate this field onto the destination grid
      if ( associated(src_data) ) deallocate(src_data)
      allocate( src_data(nx*ny) )
      src_data(1:nx*ny) = wrk(1:nx*ny)
      call regrid_Field(src_vinfo%name, src_data, dst_vinfo%name, dst_data, &
                        method=trim(regrid_method), norm_type=this_norm_opt, &
                        use_src_mask=use_src_mask, use_dst_mask=use_dst_mask, &
                        src_grid_id=id_src_grid, dst_grid_id=id_dst_grid)

      if (verbose > -1) then
        write(6,'(a,2g14.6,a,i8)')'source min,max: ', &
            minval(src_data, mask=src_data.ne.fillvalue), &
            maxval(src_data, mask=src_data.ne.fillvalue), &
            '  fillvalue count = ',count(src_data.eq.fillvalue)
        write(6,'(a,2g14.6,a,i8)')'  dest min,max: ', &
            minval(dst_data, mask=dst_data.ne.fillvalue), &
            maxval(dst_data, mask=dst_data.ne.fillvalue), &
            '  fillvalue count = ',count(dst_data.eq.fillvalue)
      endif

#undef test_remap2d
#ifdef test_remap2d
      !--- Interpolate onto the output grid using the remap2d interface

      !--- Allocate and assign the input array for remap2d
      if ( allocated(fld_in) ) deallocate( fld_in )
      allocate( fld_in(nx,ny) )
      fld_in(:,:) = reshape( wrk, (/ nx, ny /) )

      !--- Allocate the output array used in remap2d
      if ( allocated(fld_out) ) deallocate( fld_out )
      allocate( fld_out(nx2,ny2) )
      fld_out = 0.0
                     
      if (verbose > -1) then
        write(6,'(a,2g14.6,a,i8)')'source min,max: ', &
            minval(fld_in,mask=fld_in.ne.fillvalue), &
            maxval(fld_in,mask=fld_in.ne.fillvalue), &
            '  fillvalue count = ',count(fld_in.eq.fillvalue)
        write(6,'(a,2g14.6,a,i8)')'  dest min,max: ', &
            minval(fld_out,mask=fld_out.ne.fillvalue), &
            maxval(fld_out,mask=fld_out.ne.fillvalue), &
            '  fillvalue count = ',count(fld_out.eq.fillvalue)
      endif

      !--- Compare with the previous interpolation
      if ( allocated(dst_arr2) ) deallocate(dst_arr2)
      allocate( dst_arr2(nx2*ny2) )
      dst_arr2 = reshape( fld_out, (/nx2*ny2/) )
      dst_arr2 = dst_arr - dst_arr2
      write(6,'(a,2g14.6,a,i8)')'dst_arr2 min,max: ', &
          minval(dst_arr2,mask=dst_arr2.ne.fillvalue), &
          maxval(dst_arr2,mask=dst_arr2.ne.fillvalue), &
          '  fillvalue count = ',count(dst_arr2.eq.fillvalue)
      deallocate( dst_arr2 )
      deallocate( fld_in   )
      deallocate( fld_out  )
#endif

      if ( write_orig_fld ) then
        !--- Write the original input record to the output file prior
        !--- to writing the interpolated field to the same output file
        call putfld2(iuorig, wrk, ibuf, maxx)
      endif

      !--- If the interoplated field was overwritten in new_map then reset ibuf3
      select case ( trim(adjustl(overwrite_with)) )
        case ("grid_frac")
          !--- Over write the output array with the grid fraction array.
          !--- For a conservative remapping this is the area fraction of each grid
          !--- cell that participates in the remapping.
          !--- For other remapping options, the grid frac array is one where the
          !--- grid point participates in the remapping and zero otherwise.
          !--- In all cases this fraction is determined by the source grid mask, not
          !--- the grid mask of the destination array.
          ibuf(3) = transfer("FRAC",1)
        case ("grid_area")
          !--- Over write the output array with the grid area array
          !--- The grid area array is only computed during a conservative remapping
          !--- All other remapping types will fill this array with zeros
          ibuf(3) = transfer("AREA",1)
      end select

      if ( nc_out ) then
        !--- Write the interpolated field to a netcdf file
        curr_time = ibuf(2)
        nc_file = " "
        nc_file = trim(curr_file)//".nc"
        nc_name = " "
        nc_name = transfer(ibuf(3),"XXXX")
        nwrds = nx2*ny2
        if ( nrec == 1 ) then
          call write_nc1(nc_file, dst_data, nc_name, curr_time, dst_grid, nwrds, create=.true.)
        else
          call write_nc1(nc_file, dst_data, nc_name, curr_time, dst_grid, nwrds)
        endif

      else
        !--- Dump interpolated field to a CCCma binary file
        nwrds = nx2*ny2
        wrk(1:nwrds) = dst_data(1:nwrds)
        ibuf(5) = nx2
        ibuf(6) = ny2
        ibuf(8) = 1
        call putfld2(iuout, wrk, ibuf, maxx)
      endif

      if ( write_reinterp_diff ) then
        !--- Interpolate back to the original grid and compare with input data
        if ( allocated(src_arr2) ) deallocate( src_arr2 )
        allocate( src_arr2(nx*ny) )
        !--- Always ask for the output array to be returned with the cyclic
        !--- longitude included so that this twice interpolated array may be
        !--- compared with the original field that also contains a cyclic lon
        src_dims(1:2) = (/nx2,ny2/)
        if (verbose > -1) then
          write(6,*)'============= Reverse interpolation'
          write(6,'(a,2g14.6,a,i8)')'source min,max: ', &
              minval(dst_arr, mask=dst_arr.ne.fillvalue), &
              maxval(dst_arr, mask=dst_arr.ne.fillvalue), &
              '  fillvalue count = ',count(dst_arr.eq.fillvalue)
          write(6,'(a,2g14.6,a,i8)')'  dest min,max: ', &
              minval(src_arr2, mask=src_arr2.ne.fillvalue), &
              maxval(src_arr2, mask=src_arr2.ne.fillvalue), &
              '  fillvalue count = ',count(src_arr2.eq.fillvalue)
        endif
        wrk(:) = 0.0
        where (src_arr == fillvalue .or. src_arr2 == fillvalue)
          wrk(1:nx*ny) = fillvalue
        elsewhere
          wrk(1:nx*ny) = src_arr - src_arr2
        end where

        if ( nc_out ) then
          !--- Write the interpolated field to a netcdf file
          curr_time = ibuf(2)
          nc_file = " "
          nc_file = trim(curr_file)//".nc"
          nc_name = " "
          nc_name = transfer(ibuf(3),"XXXX")//"_diff"
          nwrds = nx*ny
          src_arr2(1:nx*ny) = wrk(1:nx*ny)
          call write_nc1(nc_file, src_arr2, nc_name, curr_time, src_grid, nwrds)
        else
          !--- Dump the difference between the original data and twice interpolated data
          !--- (ie interpolated onto output grid then back to input grid)
          ibuf(3) = nc4to8("DIFF")
          ibuf(5) = nx
          ibuf(6) = ny
          ibuf(8) = 1
          call putfld2(iuout, wrk, ibuf, maxx)
        endif

        deallocate( src_arr2 )
      endif

      !--- Clean up before reading the next record
      if ( allocated(src_arr) ) deallocate( src_arr )

      !--- Do not exceed the maximum number of records
      !--- that may be read from a single input file
      if (nrec .ge. max_rec_read) exit

    enddo READ_RECORDS

    !--- Close the current input file
    call cc_close(iuin,trim(curr_file),remote,1)

  enddo READ_FILES

  if ( nc_out ) then
    write(6,'(2a)')'Created file ',trim(nc_file)
  else
    close(iuout)
    write(6,*)"Create output file containing interpolated data ",trim(fout)
  endif

  if ( write_orig_fld ) then
    close(iuorig)
    write(6,*)"Create output file containing original data ",trim(forig)
  endif

  !--- Remove any links to the grid definition files
  if ( remote_src_grid ) then
    strng=' '
    strng="release "//trim(src_grid)
    call system(trim(strng))
  endif
  if ( remote_dst_grid ) then
    strng=' '
    strng="release "//trim(dst_grid)
    call system(trim(strng))
  endif

  !--- Clean up
  if ( allocated(src_arr) ) deallocate( src_arr )
  if ( allocated(dst_arr) ) deallocate( dst_arr )
  if ( associated(src_mask_values) ) deallocate(src_mask_values)
  if ( associated(dst_mask_values) ) deallocate(dst_mask_values)

  if (use_timers) call timer_cpu( 'stop', idx=1, tag='remap' )
  if (use_timers) call timer_cpu( 'print' )

  !--- Finalize ESMF
  call ESMF_Finalize()

contains

  !--- Start/Stop a given timer or print times
  subroutine timer_cpu(action, idx, tag)
    implicit none

    character(*), intent(in) :: action
    integer, intent(in), optional :: idx
    character(*), intent(in), optional :: tag

    !--- Local
    integer, parameter :: max_timers=100
    logical, save, dimension(max_timers) :: in_use = .false.
    real, save, dimension(max_timers) :: total_cpu = 0.0
    real, save, dimension(max_timers) :: time_in  = 0.0
    real, save, dimension(max_timers) :: time_out = 0.0
    character(32), save, dimension(max_timers) :: tags = " "
    real :: curr_cpu = 0.0
    integer :: curr_timer

    if ( len_trim(action) .le. 0 ) then
      write(6,'(a)')'timer_cpu: Missing action.'
      call abort
    endif

    if ( trim(adjustl(action)) .ne. 'print' ) then
      if ( .not. present(idx) ) then
        write(6,'(a)')'timer_cpu: idx must be present unless action=print.'
        call abort
      endif

      if ( idx.le.0 .or. idx.gt.max_timers ) then
        write(6,'(a,i6,a)')'timer_cpu: Index = ',idx,' is out of range.'
        call abort
      endif

      if ( present(tag) ) then
        tags(idx) = trim(adjustl(tag))
      endif
    endif

    select case ( trim(adjustl(action)) )
      case ('start')
        if ( in_use(idx) ) then
          write(6,'(a,i4,a)')'timer_cpu: Attempting to start timer ',idx, &
            ' but it is currently running ...Ignored'
          return
        endif
        in_use(idx) = .true.
        call cpu_time( time_in(idx) )

      case ('stop', 'stop_and_print')
        if ( .not. in_use(idx) ) then
          write(6,'(a,i4,a)')'timer_cpu: Attempting to stop timer ',idx, &
            ' but it is not running ...Ignored'
          return
        endif
        call cpu_time( time_out(idx) )
        curr_cpu = time_out(idx) - time_in(idx)
        total_cpu(idx) = total_cpu(idx) + curr_cpu
        if ( trim(adjustl(action)) .eq. 'stop_and_print' ) then
          if ( present(tag) ) then
            write(6,'(a,i4,a,g16.8,a,g16.8,4x,a)') 'Stopped timer ',idx,' after ', &
              curr_cpu,' seconds.  Total=',total_cpu(idx),trim(tag)
          else
            write(6,'(a,i4,a,g16.8,a,g16.8)') 'Stopped timer ',idx,' after ', &
              curr_cpu,' seconds.  Total=',total_cpu(idx)
          endif
        endif
        in_use(idx) = .false.

      case ('print')
        do curr_timer=1,max_timers
          if ( total_cpu(curr_timer) .gt. 0.0 ) then
            if ( len_trim(tags(curr_timer)).gt.0 ) then
              write(6,'(a,i4,a,g16.8,a,4x,a)')'Total for timer ',curr_timer, &
                ' is ',total_cpu(curr_timer),' seconds.',trim(tags(curr_timer))
            else
              write(6,'(a,i4,a,g16.8,a)')'Total for timer ',curr_timer, &
                ' is ',total_cpu(curr_timer),' seconds.'
            endif
          endif
        enddo

      case default
        write(6,*)"timer_cpu: Invalid action ",trim(action)
        call abort
    end select

  end subroutine timer_cpu

  subroutine def_input(iopt, fin, lat, lon, nx, ny, overlap, imask)

    implicit none

    real :: fin(nx,ny)
    real(8) :: lat(:), lon(:)
    integer :: iopt, nx, ny
    logical :: overlap
    integer, optional :: imask(:)

    !--- local
    integer :: nxx, ix, iy, tic, toc
    integer, allocatable :: imask2d(:,:)

    !xxx write(6,*)'lat:'
    !xxx write(6,'(10f10.4)')lat(1:ny)
    !xxx write(6,*)'lon:'
    !xxx write(6,'(10f10.4)')lon(1:nx)

    if (overlap) then
      nxx = nx-1
    else
      nxx = nx
    endif

    !--- imask will always be defined on a non-overlapping grid
    allocate( imask2d(nxx,ny) )

    if ( present(imask) ) then
      !--- assign imask to a local 2d array
      imask2d = reshape( imask, (/ nxx, ny /) )
    else
      imask2d = 1
    endif

    select case (iopt)
      case (1)
        !--- Create a checker board pattern of 1's and 3's
        !--- Apply the input integer mask to these values
        tic=1
        do iy=1,ny
          tic=-tic
          toc=tic
          do ix=1,nxx
            toc=-toc
            if (toc.gt.0) then
              fin(ix,iy) = 1.0 * imask2d(ix,iy)
            else
              fin(ix,iy) = 3.0 * imask2d(ix,iy)
            endif
          enddo
        enddo

      case (2)
        !--- pseudo-spherical harmonic l=2,m=2
        do iy=1,ny
          do ix=1,nxx
            fin(ix,iy) = imask2d(ix,iy) * &
              (2.d0 + cos(lat(iy))**2 * cos(2.d0*lon(ix)))
          enddo
        enddo

      case default
        write(6,'(2a,i3)')'def_input: Invalid iopt = ',iopt
        call xit("def_input",-1)
    end select

    if (overlap) then
      do iy=1,ny
        fin(nx,iy) = fin(1,iy)
      enddo
    endif

    deallocate( imask2d )

  end subroutine def_input

  subroutine read_grid_dims(nx, ny, grid_file, lat, lon, imask, msg)

    use netcdf

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

    !--- Name of file containing the grid info
    character(*), intent(in) :: grid_file

    integer, intent(out), optional :: imask(:)
    real(8), intent(out), optional :: lat(:)
    real(8), intent(out), optional :: lon(:)
    character(*), intent(out), optional :: msg

    !--- local

    character(128) :: title
    character(256) :: etxt

    integer (4) :: ncstat, ncfid, ncvid, grid_rank, grid_size
    integer (4), allocatable, dimension(:) :: grid_dims
    real(8), allocatable :: grid_lon(:,:), grid_lat(:,:), wrk(:)

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

    !--- Open the grid info file
    ncstat = nf90_open(grid_file, NF90_NOWRITE, ncfid)
    etxt = " "
    write(etxt,'("read_grid_dims: Problem opening file ",a)')trim(grid_file)
    call nc_error_handler(ncfid, ncstat, 1, msg=trim(etxt))

    !--- Read and assign the "title" attribute
    title = ' '
    ncstat = nf90_get_att(ncfid,NF90_GLOBAL,'title',title)
    call nc_error_handler(ncfid, ncstat, 2)

    !--- Read the size of this grid
    ncstat = nf90_inq_dimid(ncfid, 'grid_size', ncvid)
    call nc_error_handler(ncfid, ncstat, 3)
    ncstat = nf90_inquire_dimension(ncfid, ncvid, len=grid_size)
    call nc_error_handler(ncfid, ncstat, 4)

    !--- Read number of dims in this grid
    ncstat = nf90_inq_dimid(ncfid, 'grid_rank', ncvid)
    call nc_error_handler(ncfid, ncstat, 3)
    ncstat = nf90_inquire_dimension(ncfid, ncvid, len=grid_rank)
    call nc_error_handler(ncfid, ncstat, 4)

    allocate( grid_dims( grid_rank ) )

    !--- Read the shape of the source grid
    ncstat = nf90_inq_varid(ncfid, 'grid_dims', ncvid)
    call nc_error_handler(ncfid, ncstat, 5)
    ncstat = nf90_get_var(ncfid, ncvid, grid_dims)
    call nc_error_handler(ncfid, ncstat, 6)

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

    deallocate( grid_dims )

    if ( present(lat) ) then
      !--- Also read grid center latitude
      ncstat = nf90_inq_varid(ncfid, 'grid_center_lat', ncvid)
      call nc_error_handler(ncfid, ncstat, 9)
      !--- Read lat as an array of shape (nx,ny)
      allocate( grid_lat(nx,ny), wrk(nx*ny) )
      ncstat = nf90_get_var(ncfid, ncvid, wrk)
      call nc_error_handler(ncfid, ncstat, 10)
      !--- Return lat as a 1D array of length ny
      grid_lat = reshape( wrk, (/nx,ny/) )
      lat(1:ny) = grid_lat(1,1:ny)
      deallocate( grid_lat, wrk )
    endif

    if ( present(lon) ) then
      !--- Also read grid center longitude
      ncstat = nf90_inq_varid(ncfid, 'grid_center_lon', ncvid)
      call nc_error_handler(ncfid, ncstat, 11)
      !--- Read lon as an array of shape (nx,ny)
      allocate( grid_lon(nx,ny), wrk(nx*ny) )
      ncstat = nf90_get_var(ncfid, ncvid, wrk)
      call nc_error_handler(ncfid, ncstat, 12)
      !--- Return lon as a 1D array of length nx
      grid_lon = reshape( wrk, (/nx,ny/) )
      lon(1:nx) = grid_lon(1:nx,1)
      deallocate( grid_lon, wrk )
    endif

    if ( present(imask) ) then
      !--- Also read the integer mask
      ncstat = nf90_inq_varid(ncfid, 'grid_imask', ncvid)
      call nc_error_handler(ncfid, ncstat, 7)
      ncstat = nf90_get_var(ncfid, ncvid, imask)
      call nc_error_handler(ncfid, ncstat, 8)
    endif

    !--- Close the netcdf file
    ncstat = nf90_close(ncfid)
    call nc_error_handler(ncfid, ncstat, 13)

    if ( present(msg) ) then
      if ( trim(adjustl(msg)) .eq. "write" ) then
        !--- Write a diagnostic message to stdout
        write(6,'(a,2i6,2a)')' shape: ',nx,ny,'   title: ',trim(title)
      else
        !--- Write a diagnostic message into a string to pass back to caller
        write(msg,'(a,2i6,2a)')' shape: ',nx,ny,'   title: ',trim(title)
      endif
    endif

  end subroutine read_grid_dims

  subroutine write_nc1(fname, var, vname, curr_time, grid, nwrds, create)

    use netcdf

    !--- Name of file to create
    character(*), intent(in) :: fname

    !--- Variable values in a 1D array
    real(kind=8), intent(in) :: var(:)

    !--- A string to use as the variable name in the netcdf output file
    character(*), intent(in) :: vname

    !--- The time associated with this variable
    integer(kind=8), intent(in) :: curr_time

    !--- Name of file containing grid information
    character(*), intent(in) :: grid

    !--- The number of words in the data passed in var
    integer, intent(in) :: nwrds

    !--- Name of file containing grid information
    logical, intent(in), optional :: create

    !--- local

    !--- Latitude and longitude values
    real(kind=8), allocatable :: lon(:)
    real(kind=8), allocatable :: lat(:)
    real(kind=8), allocatable :: time_vals(:)

    logical :: create_new_file, found, new_lat, new_lon
    real(kind=8), allocatable :: var_tmp(:)
    character(len=NF90_MAX_NAME) :: curr_name
    character(len=2) :: ch2
    integer :: nx_grid, ny_grid, lon_n, lat_n, dim_n
    integer :: ix, iy, nxx, idx1, idx2
    integer (kind=4) :: ncstat, fid, vid, idx
    integer (kind=4) :: ndims, nvars, natts, recdim_id, dimsz
    integer (kind=4), save :: nx, ny, nt
    integer (kind=4), save :: lon_dim_id, lat_dim_id, time_dim_id
    integer (kind=4), save :: lon_id, lat_id, time_id

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

    if (len_trim(fname) <= 0) then
      write(6,'(a)')'write_nc1: fname is not defined.'
      call flush(6)
      call xit("write_nc1",-1)
    endif

    if (len_trim(vname) <= 0) then
      write(6,'(a)')'write_nc1: vname is not defined.'
      call flush(6)
      call xit("write_nc1",-2)
    endif

    !--- Determine lat/lon cooordinates from the grid description file
    call read_grid_dims(nx_grid, ny_grid, grid)
    allocate( lon(nx_grid), lat(ny_grid) )
    call read_grid_dims(nx_grid, ny_grid, grid, lon=lon, lat=lat)

    nx = nx_grid
    ny = ny_grid

    if ( nx<=0 .or. ny<=0 ) then
      write(6,'(a)')'write_nc1: Invalid dimension sizes.'
      write(6,'(a,2i8)')'  nx, ny = ',nx,ny
      call flush(6)
      call xit("write_nc1",-3)
    endif

    if ( present(create) ) then
      !--- Use the input value of create
      create_new_file = create
    else
      !--- If the file exists then read it, if not, create a new file
      inquire(file=trim(adjustl(fname)), exist=found)
      if (found) then
        create_new_file = .false.
      else
        create_new_file = .true.
      endif
    endif

    if ( create_new_file ) then

      !--- Create the output netcdf file, overwriting any file by that name
      ncstat = nf90_create(fname, NF90_CLOBBER, fid)
      call nc_error_handler(fid, ncstat, 6)

      !--- Define dimension IDs
      ncstat = nf90_def_dim(fid, "time", NF90_UNLIMITED, time_dim_id)
      call nc_error_handler(fid, ncstat, 6)

      ncstat = nf90_def_dim(fid, "lon1", nx, lon_dim_id)
      call nc_error_handler(fid, ncstat, 4)

      ncstat = nf90_def_dim(fid, "lat1", ny, lat_dim_id)
      call nc_error_handler(fid, ncstat, 5)

      !--- Define dimension variables
      ncstat = nf90_def_var(fid, "time", NF90_DOUBLE, time_dim_id, time_id)
      call nc_error_handler(fid, ncstat, 9)

      ncstat = nf90_put_att(fid, time_id, "long_name", "Time")
      call nc_error_handler(fid, ncstat, 7)

      ncstat = nf90_put_att(fid, time_id, "units", "ibuf2")
      call nc_error_handler(fid, ncstat, 7)

      ncstat = nf90_def_var(fid, "lon1", NF90_DOUBLE, lon_dim_id, lon_id)
      call nc_error_handler(fid, ncstat, 7)

      ncstat = nf90_put_att(fid, lon_id, "long_name", "Longitude")
      call nc_error_handler(fid, ncstat, 7)

      ncstat = nf90_put_att(fid, lon_id, "units", "degrees_east")
      call nc_error_handler(fid, ncstat, 7)

      ncstat = nf90_def_var(fid, "lat1", NF90_DOUBLE, lat_dim_id, lat_id)
      call nc_error_handler(fid, ncstat, 8)

      ncstat = nf90_put_att(fid, lat_id, "long_name", "Latitude")
      call nc_error_handler(fid, ncstat, 7)

      ncstat = nf90_put_att(fid, lat_id, "units", "degrees_north")
      call nc_error_handler(fid, ncstat, 7)

      !--- Define the output variable and its attributes
      ncstat = nf90_def_var(fid, trim(adjustl(vname)), NF90_FLOAT, &
                            (/lon_dim_id, lat_dim_id, time_dim_id/), vid)
      call nc_error_handler(fid, ncstat, 10)

      !--- Get out of define mode
      ncstat = nf90_enddef(fid)
      call nc_error_handler(fid, ncstat, 11)

      !--- Write coordinate variable data
      ncstat = nf90_put_var(fid, lon_id, lon(1:nx))
      call nc_error_handler(fid, ncstat, 12)

      ncstat = nf90_put_var(fid, lat_id, lat(1:ny))
      call nc_error_handler(fid, ncstat, 13)

      nt = 0

    else
      !--- Open an existing netcdf file for writing in data mode
      ncstat = nf90_open(fname, NF90_WRITE, fid)
      call nc_error_handler(fid, ncstat, 11)

      !--- Find out about existing variables
      ncstat = nf90_Inquire(fid, ndims, nvars, natts, time_dim_id)
      call nc_error_handler(fid, ncstat, 11)

      !--- Determine the current size of the record dimension
      ncstat = nf90_Inquire_Dimension(fid, time_dim_id, name=curr_name, len=dimsz)
      call nc_error_handler(fid, ncstat, 11)
      nt = dimsz

      !--- Determine the variable ID of the time coordinate
      ncstat = nf90_inq_varid(fid, "time", time_id)
      call nc_error_handler(fid, ncstat, 11)

      !--- Determine if curr_time has already been written and if so what is nt
      if ( allocated(time_vals) ) deallocate(time_vals)
      allocate( time_vals(dimsz) )
      ncstat = nf90_get_var(fid, time_id, time_vals)
      call nc_error_handler(fid, ncstat, 11)
      do idx = 1,dimsz
        if ( abs(curr_time-time_vals(idx)) .lt. &
             5.0*abs(nearest(time_vals(idx),1.0) - time_vals(idx)) ) then
          !--- We have a match
          !--- Note nt is incremented below prior to writing
          nt = idx - 1
          exit
        endif
      enddo
      deallocate(time_vals)

      !--- Put the file into define mode
      ncstat = nf90_redef(fid)
      call nc_error_handler(fid, ncstat, 10)

      !--- If a dimension is named lon?? or lat?? and is the same size as nx or ny
      !--- then assume that it is the dimension to be used for this variable
      lon_dim_id = -1
      lat_dim_id = -1
      lon_n = 0
      lat_n = 0
      do idx=1,ndims
        ncstat = nf90_Inquire_Dimension(fid, idx, name=curr_name, len=dimsz)
        call nc_error_handler(fid, ncstat, 11)
        !--- Ensure any NULL in curr_name(4:5) is converted to a space
        if ( curr_name(4:4) .eq. achar(0) ) curr_name(4:4) = " "
        if ( curr_name(5:5) .eq. achar(0) ) curr_name(5:5) = " "
        if ( curr_name(1:3) .eq. "lon" ) then
          ch2 = curr_name(4:5)
          read(ch2,*)dim_n
          lon_n = max(lon_n, dim_n)
          if ( dimsz == nx ) lon_dim_id = idx
        endif
        if ( curr_name(1:3) .eq. "lat" ) then
          ch2 = curr_name(4:5)
          read(ch2,*)dim_n
          lat_n = max(lat_n, dim_n)
          if ( dimsz == ny ) lat_dim_id = idx
        endif
      enddo

      new_lon = .false.
      if ( lon_dim_id < 0 ) then
        !--- The input lon dimension was not found in the file
        new_lon = .true.
        if ( lon_n == 0 ) then
          lon_n = 1
        else
          lon_n = lon_n + 1
        endif
        ch2 = " "
        write(ch2,'(i2)')lon_n
        curr_name = " "
        curr_name = "lon"//trim(adjustl(ch2))
        !--- Define the new lon dimension ID
        ncstat = nf90_def_dim(fid, trim(curr_name), nx, lon_dim_id)
        call nc_error_handler(fid, ncstat, 4)

        !--- Define the new lon coordinate variable
        ncstat = nf90_def_var(fid, trim(curr_name), NF90_DOUBLE, lon_dim_id, lon_id)
        call nc_error_handler(fid, ncstat, 7)

        ncstat = nf90_put_att(fid, lon_id, "long_name", "Longitude")
        call nc_error_handler(fid, ncstat, 7)

        ncstat = nf90_put_att(fid, lon_id, "units", "degrees_east")
        call nc_error_handler(fid, ncstat, 7)
      endif

      new_lat = .false.
      if ( lat_dim_id < 0 ) then
        !--- The input lat dimension was not found in the file
        new_lat = .true.
        if ( lat_n == 0 ) then
          lat_n = 1
        else
          lat_n = lat_n + 1
        endif
        ch2 = " "
        write(ch2,'(i2)')lat_n
        curr_name = " "
        curr_name = "lat"//trim(adjustl(ch2))
        !--- Define the new lat dimension ID
        ncstat = nf90_def_dim(fid, trim(curr_name), ny, lat_dim_id)
        call nc_error_handler(fid, ncstat, 4)

        !--- Define the new lat coordinate variable
        ncstat = nf90_def_var(fid, trim(curr_name), NF90_DOUBLE, lat_dim_id, lat_id)
        call nc_error_handler(fid, ncstat, 7)

        ncstat = nf90_put_att(fid, lat_id, "long_name", "Latitude")
        call nc_error_handler(fid, ncstat, 7)

        ncstat = nf90_put_att(fid, lat_id, "units", "degrees_east")
        call nc_error_handler(fid, ncstat, 7)
      endif

      vid = -1
      do idx=1,nvars
        ncstat = nf90_Inquire_Variable(fid, idx, name=curr_name)
        call nc_error_handler(fid, ncstat, 11)
        if ( trim(adjustl(curr_name)) .eq. trim(adjustl(vname)) ) then
          !--- The input variable already exists in this file
          vid = idx
          exit
        endif
      enddo

      if ( vid<0 ) then
        !--- The input variable was not found in the file
        ncstat = nf90_def_var(fid, trim(adjustl(vname)), NF90_FLOAT, &
                              (/lon_dim_id, lat_dim_id, time_dim_id/), vid)
        call nc_error_handler(fid, ncstat, 10)
      endif

      !--- Get out of define mode
      ncstat = nf90_enddef(fid)
      call nc_error_handler(fid, ncstat, 11)

      if ( new_lon ) then
        !--- Write lon coordinate variable data if newly defined
        ncstat = nf90_put_var(fid, lon_id, lon(1:nx))
        call nc_error_handler(fid, ncstat, 12)
      endif

      if ( new_lat ) then
        !--- Write lat coordinate variable data if newly defined
        ncstat = nf90_put_var(fid, lat_id, lat(1:ny))
        call nc_error_handler(fid, ncstat, 13)
      endif

    endif

    !--- Increment the time counter
    nt = nt + 1

    ncstat = nf90_put_var( fid, time_id, curr_time, start=(/nt/) )
    call nc_error_handler(fid, ncstat, 14)

    if ( allocated(var_tmp) ) deallocate(var_tmp)
    allocate( var_tmp(nx*ny) )
    if ( nwrds == nx*ny ) then
      !---- The input data does not contain an overlap lon
      var_tmp(1:nx*ny) = var(1:nx*ny)
    else if (nwrds == (nx+1)*ny) then
      !---- The input data contains an overlap lon
      nxx = nx + 1
      idx2=0
      do iy=1,ny
        do ix=1,nxx
          !--- idx1 is the index into the overlap array
          idx1 = (iy-1)*nxx + ix
          if (mod(idx1,nxx).eq.0) then
            !--- This is the overlap longitude ...ignore it
            cycle
          endif
          idx2=idx2+1
          var_tmp(idx2) = var(idx1)
        end do
      end do
      if (idx2 .ne. nx*ny) then
        write(6,'(a)')'write_nc1: Problem removing cyclic longitude from input arr.'
        write(6,'(a,i8)')'  Grid size after removal =  ',idx2
        write(6,'(a,i8)')'       Grid size expected =  ',nx*ny
        call xit("write_nc1",-90)
      endif
    else
      write(6,'(a)')'write_nc1: Invalid input array size.'
      write(6,'(a,2i8)')'  nx, ny = ',nx,ny
      write(6,'(a,i8)')'   nwrds = ',nwrds
      call flush(6)
      call xit("write_nc1",-91)
    endif
    ncstat = nf90_put_var( fid, vid, var_tmp, start=(/1_4,1_4,nt/), count=(/nx,ny,1_4/) )
    call nc_error_handler(fid, ncstat, 92)

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

    !--- Clean up
    if ( allocated(lon)     ) deallocate(lon)
    if ( allocated(lat)     ) deallocate(lat)
    if ( allocated(var_tmp) ) deallocate(var_tmp)

  end subroutine write_nc1

end program regrid
