ESMF_FieldRegridXGUTest Program

      program ESMF_FieldRegridXGUTest

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

#include "ESMF.h"

!==============================================================================
!BOPI
! !PROGRAM: ESMF_FieldRegridXGUTest - Unit tests for Field RegridXG Testing
!
! !DESCRIPTION:
!
! The code in this file drives F90 Field RegridXG Testing unit tests.
! The companion folder Field\/src contains the definitions for the
! Field methods.
!EOPI
!-----------------------------------------------------------------------------
! !USES:
    use ESMF_TestMod     ! test methods
    use ESMF

    implicit none

!------------------------------------------------------------------------------
! The following line turns the CVS identifier string into a printable variable.
    character(*), parameter :: version = &
      '$Id$'

    ! cumulative result: count failures; no failures equals "all pass"
    integer :: result = 0

    ! individual test result code
    integer :: rc = 1

    ! individual test failure message
    character(ESMF_MAXSTR) :: failMsg
    character(512) :: name

    call ESMF_TestStart(ESMF_SRCLINE, rc=rc)
    if (rc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT)
 
#ifdef ESMF_TESTEXHAUSTIVE

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xg_online(3,3,4,4,1.,1.,0.6,0.6,tag='small test', &
      maxnpet=2, indexflag=ESMF_INDEX_GLOBAL, rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid online and regrid through xgrid, overlapping cut"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xg_online(4,4,3,3,0.6,0.6,1.,1.,tag='reverse small test', &
      maxnpet=2, indexflag=ESMF_INDEX_GLOBAL, rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid online and regrid through xgrid, overlapping cut"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xg_online(10,10,14,14,0.1,0.1,0.06,0.06,tag='medium size test A1', &
      maxnpet=4, indexflag=ESMF_INDEX_GLOBAL, rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid online and regrid through xgrid, overlapping cut"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xg_online(10,10,10,10,0.1,0.1,0.06,0.06,tag='medium size test A2', &
      maxnpet=4, indexflag=ESMF_INDEX_GLOBAL, rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid online and regrid through xgrid, overlapping cut"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xg_online(10,10,10,10,0.1,0.1,0.1,0.1,tag='medium size test A3', &
      maxnpet=4, indexflag=ESMF_INDEX_GLOBAL, rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid online and regrid through xgrid, overlapping cut"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xg_online(10,10,14,14,0.1,0.1,0.06,0.06,tag='reverse medium size test', &
      maxnpet=4, indexflag=ESMF_INDEX_GLOBAL, rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid online and regrid through xgrid, overlapping cut"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xg_online(3,4,4,6,1.,1.,0.6,0.6,tag='uneven small test', &
      maxnpet=2, indexflag=ESMF_INDEX_GLOBAL, rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid online and regrid through xgrid, overlapping cut"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xg_online(4,6,3,4,0.6,0.6,1.,1.,tag='reverse uneven small test', &
      maxnpet=2, indexflag=ESMF_INDEX_GLOBAL, rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid online and regrid through xgrid, overlapping cut"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xg_online(80,80,40,40,0.5,0.5,1.,1.,tag='exact large test A2', &
      maxnpet=8, indexflag=ESMF_INDEX_GLOBAL, rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid online and regrid through xgrid, overlapping cut"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xg_online(40,40,80,80,1.,1.,0.5,0.5,tag='reverse exact large test', &
      maxnpet=8, indexflag=ESMF_INDEX_GLOBAL, rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid online and regrid through xgrid, overlapping cut"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xgSph(6,5,6,4,1.,1.,1.,1.,165.,30.,165.,30., &
      maxnpet=2, tag='small regional sphere exact', rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid and regrid through xgrid, spherical grids"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xgSph(6,5,6,4,1.,1.,1.,1.,-165.,30.,-165.,30., &
      maxnpet=2, tag='small regional sphere exact negative lon', rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid and regrid through xgrid, spherical grids"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xgSph(5,6,5,6,2.,1.,1.,1.,-165.,30.,-165.,30., &
      maxnpet=2, tag='small regional sphere contain', rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid and regrid through xgrid, spherical grids"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xgSph(10,20,20,20,1.,1.,0.5,1.,-165.,30.,-165.,30., &
      maxnpet=8, tag='medium regional sphere cut B1', rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid and regrid through xgrid, spherical grids"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xgSph(10,20,20,20,1.,1.,0.6,1.,-165.,30.,-165.,30., &
      maxnpet=8, tag='medium regional sphere overlap', rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid and regrid through xgrid, spherical grids"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xgSph(10,20,20,20,1.,1.,0.6,1.,-165.,30.,-166.,30., &
      maxnpet=8, tag='medium regional sphere overlap2', rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid and regrid through xgrid, spherical grids"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xgSph(40,40,80,80,1.,1.,0.6,0.6,-165.,30.,-168.,25., &
      tag='large regional sphere overlap B2', rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid and regrid through xgrid, spherical grids"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xgSph(80,80,40,40,0.6,0.6,1.,1.,-168.,25.,-165.,30., &
      tag='reverse large regional sphere overlap', rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid and regrid through xgrid, spherical grids"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)


    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xgSph(45,45,90,90,8.,4.,4.,2.,-180.,-90.,0.,-90., &
      scheme=ESMF_REGRID_SCHEME_FULL3D, tag='large full sphere cut C1', rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid and regrid through xgrid, full spherical grids"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xgSph(90,90,45,45,4.,2.,8.,4.,0.,-90.,-180.,-90., &
      scheme=ESMF_REGRID_SCHEME_FULL3D, &
      tag='reverse large full sphere cut', rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid and regrid through xgrid, full spherical grids"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xgSph(60,60,90,90,6.,3.,4.,2.,-180.,-90.,0.,-90., &
      scheme=ESMF_REGRID_SCHEME_FULL3D, tag='large full sphere latlonclip C2', rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid and regrid through xgrid, full spherical grids"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)

    !------------------------------------------------------------------------
    !EX_UTest
    call test_regrid2xgSph(90,90,60,60,4.,2.,6.,3.,0.,-90.,-180.,-90., &
      scheme=ESMF_REGRID_SCHEME_FULL3D, &
      tag='reverse large full sphere latlonclip', rc=rc)
    write(failMsg, *) ""
    write(name, *) "Regrid then create xgrid and regrid through xgrid, full spherical grids"
    call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)
#endif
    call ESMF_TestEnd(ESMF_SRCLINE)

contains 
#define ESMF_METHOD "ESMF_TESTS"

!------------------------------------------------------------------------
  subroutine test_regrid2xg_online(atm_nx, atm_ny, ocn_nx, ocn_ny, atm_dx, atm_dy, &
    ocn_dx, ocn_dy, atm_sx, atm_sy, ocn_sx, ocn_sy, tag, maxnpet, minnpet, indexflag, rc)
    ! arguments
    integer, intent(in)                       :: atm_nx, atm_ny, ocn_nx, ocn_ny
    real(ESMF_KIND_R4), intent(in)            :: atm_dx, atm_dy, ocn_dx, ocn_dy
    real(ESMF_KIND_R4), intent(in), optional  :: atm_sx, atm_sy, ocn_sx, ocn_sy
    character(len=*), intent(in), optional    :: tag
    integer, intent(in) , optional            :: maxnpet, minnpet
    type(ESMF_INDEX_FLAG), intent(in), optional :: indexflag
    integer, intent(out), optional            :: rc

    ! local variables
    type(ESMF_Grid)                 :: grid_atm, grid_ocn
    type(ESMF_Field)                :: f_atm, f_ocn, f_xgrid, f_tmp
    real(ESMF_KIND_R8)              :: startx, starty
    integer                         :: localrc, npet, i, j, lpet
    real(ESMF_KIND_R8), pointer     :: weights(:)
    integer(ESMF_KIND_I4), pointer  :: indices(:,:)
    real(ESMF_KIND_R8), pointer     :: coordX(:), coordY(:)
    type(ESMF_XGrid)                :: xgrid
    type(ESMF_XGridSpec)            :: sparseMatA2X(1)
    integer                         :: gn(2), simax, dimax, l_minnpet, l_maxnpet
    real(ESMF_KIND_R8), pointer     :: atm(:,:), ocn(:,:), exf(:), xArea(:), xFrac(:)
    real(ESMF_KIND_R8), pointer     :: srcFracPtr(:,:), dstFracPtr(:,:)
    real(ESMF_KIND_R8), pointer     :: srcAreaPtr(:,:), dstAreaPtr(:,:)
    type(ESMF_RouteHandle)          :: rh_a2o, rh_o2a, rh_a2x, rh_x2o, rh_o2x, rh_x2a
    type(ESMF_DistGrid)             :: distgridM, distgrid1, distgrid2
    type(ESMF_XGridSpec)            :: sparseMat(1)

    type(ESMF_Field)                :: fa_atm, fa_ocn, fa_xgrid
    type(ESMF_Field)                :: aa_atm, aa_ocn, aa_xgrid
    type(ESMF_Field)                :: srcFrac, dstFrac, srcArea, dstArea
    type(ESMF_Mesh)                 :: mesh_atm, mesh_ocn, mesh_xgrid
    real(ESMF_KIND_R8)              :: srcsum(3), allsrcsum(3), dstFlux_reg, dstFlux, error
    integer                         :: eleCount, totCount(1)
    integer :: srcTermProcessing, pipeLineDepth    
    type(ESMF_VM)   :: vm

    localrc = ESMF_SUCCESS
    if(present(rc)) rc = ESMF_SUCCESS

    !------------------------------------
    ! Prepare for the calculation
    !------------------------------------
    call ESMF_VMGetCurrent(vm, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_VMGet(vm, petCount=npet, localPet=lpet, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    ! Don't run this test if there isn't enough pets or too many cuts
    l_minnpet = 1
    l_maxnpet = 128
    if(present(minnpet)) l_minnpet = minnpet
    if(present(maxnpet)) l_maxnpet = maxnpet
    if(npet > l_maxnpet .or. npet < l_minnpet) return

    if(lpet == 0 .and. present(tag)) then
      print *, '!------------------------------------'
      print *, '!', tag
      print *, '!------------------------------------'
    endif

    !------------- ATM ------------------
    ! atm grid, horizontally decomposed
    !------------------------------------
    grid_atm = ESMF_GridCreateNoPeriDim(maxIndex=(/atm_nx, atm_ny/), &
      coordSys=ESMF_COORDSYS_CART, &
      indexflag=indexflag, &
      coordDep1=(/1/), &
      coordDep2=(/2/), &
      rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_GridAddCoord(grid_atm, staggerloc=ESMF_STAGGERLOC_CENTER, &
        rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_GridAddCoord(grid_atm, staggerloc=ESMF_STAGGERLOC_CORNER, &
        rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    ! global indexing
    ! atm grid is not decomposed in the y direction
    !startx = lpet*atm_nx/npet*atm_dx
    startx = 0.
    starty = 0.
    ! compute coord
    ! X center
    call ESMF_GridGetCoord(grid_atm, localDE=0, staggerLoc=ESMF_STAGGERLOC_CENTER, &
        coordDim=1, farrayPtr=coordX, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    do i = lbound(coordX,1), ubound(coordX,1)
      coordX(i) = startx + atm_dx/2. + (i-1)*atm_dx
    enddo
    print *, 'coordX: ', coordX
    ! X corner
    call ESMF_GridGetCoord(grid_atm, localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, &
        coordDim=1, farrayPtr=coordX, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    do i = lbound(coordX,1), ubound(coordX,1)
      coordX(i) = startx + (i-1)*atm_dx
    enddo
    !print *, 'startx: ', startx, lbound(coordX, 1), 'coordX: ', coordX
    ! Y center
    call ESMF_GridGetCoord(grid_atm, localDE=0, staggerLoc=ESMF_STAGGERLOC_CENTER, &
        coordDim=2, farrayPtr=coordY, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    do i = lbound(coordY,1), ubound(coordY,1)
      coordY(i) = starty + atm_dy/2. + (i-1)*atm_dy
    enddo
    ! Y corner
    call ESMF_GridGetCoord(grid_atm, localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, &
        coordDim=2, farrayPtr=coordY, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    do i = lbound(coordY,1), ubound(coordY,1)
      coordY(i) = starty + (i-1)*atm_dy
    enddo

    !------------- OCN ------------------
    ! ocn grid, horizontally decomposed
    !------------------------------------
    grid_ocn = ESMF_GridCreateNoPeriDim(maxIndex=(/ocn_nx, ocn_ny/), &
      coordSys=ESMF_COORDSYS_CART, &
      indexflag=indexflag, &
      !gridEdgeLWidth=(/0,0/), gridEdgeUWidth=(/0,0/), &
      !regDecomp=(/npet, 1/), &
      coordDep1=(/1/), &
      coordDep2=(/2/), &
      rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_GridAddCoord(grid_ocn, staggerloc=ESMF_STAGGERLOC_CENTER, &
        rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_GridAddCoord(grid_ocn, staggerloc=ESMF_STAGGERLOC_CORNER, &
        rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    ! ocn grid is not decomposed in the y direction
    !startx = lpet*ocn_nx/npet*ocn_dx
    startx = 0.
    starty = 0.
    ! compute coord
    ! X center
    call ESMF_GridGetCoord(grid_ocn, localDE=0, staggerLoc=ESMF_STAGGERLOC_CENTER, &
        coordDim=1, farrayPtr=coordX, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    do i = lbound(coordX,1), ubound(coordX,1)
      coordX(i) = startx + ocn_dx/2. + (i-1)*ocn_dx
    enddo
    ! X corner
    call ESMF_GridGetCoord(grid_ocn, localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, &
        coordDim=1, farrayPtr=coordX, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    do i = lbound(coordX,1), ubound(coordX,1)
      coordX(i) = startx + (i-1)*ocn_dx
    enddo
    ! Y center 
    call ESMF_GridGetCoord(grid_ocn, localDE=0, staggerLoc=ESMF_STAGGERLOC_CENTER, &
        coordDim=2, farrayPtr=coordY, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    do i = lbound(coordY,1), ubound(coordY,1)
      coordY(i) = starty + ocn_dy/2. + (i-1)*ocn_dy
    enddo
    ! Y corner
    call ESMF_GridGetCoord(grid_ocn, localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, &
        coordDim=2, farrayPtr=coordY, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    do i = lbound(coordY,1), ubound(coordY,1)
      coordY(i) = starty + (i-1)*ocn_dy
    enddo

    !------------------------------------
    ! build Fields on the Grids
    !------------------------------------
    f_atm = ESMF_FieldCreate(grid_atm, typekind=ESMF_TYPEKIND_R8, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_FieldFill(f_atm, dataFillScheme="one", member=1, step=1, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    f_ocn = ESMF_FieldCreate(grid_ocn, typekind=ESMF_TYPEKIND_R8, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    srcFrac = ESMF_FieldCreate(grid_atm, typekind=ESMF_TYPEKIND_R8, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    dstFrac = ESMF_FieldCreate(grid_ocn, typekind=ESMF_TYPEKIND_R8, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    srcArea = ESMF_FieldCreate(grid_atm, typekind=ESMF_TYPEKIND_R8, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    dstArea = ESMF_FieldCreate(grid_ocn, typekind=ESMF_TYPEKIND_R8, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_FieldRegridStore(srcField=f_atm, dstField=f_ocn, &
      regridmethod=ESMF_REGRIDMETHOD_CONSERVE, &
      routehandle=rh_a2o, &
      unmappedaction = ESMF_UNMAPPEDACTION_IGNORE, &
      srcFracField=srcFrac, dstFracField=dstFrac, & 
      factorIndexList=indices, factorList=weights, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    ! make sure the numbers are consistent
    !print *, lpet, 'weights', size(weights)
    !print *, lpet, 'indices', size(indices,1),'-', size(indices,2)
    !do j = 1, size(indices,1)
    !     print *, indices(j,1), '->', indices(j,2), weights(j)
    !enddo

    call ESMF_FieldRegridStore(srcField=f_ocn, dstField=f_atm, &
      regridmethod=ESMF_REGRIDMETHOD_CONSERVE, &
      routehandle=rh_o2a, &
      unmappedaction = ESMF_UNMAPPEDACTION_IGNORE, &
      factorIndexList=indices, factorList=weights, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    ! make sure the numbers are consistent
    !print *, lpet, 'weights', size(weights)
    !print *, lpet, 'indices', size(indices,1),'-', size(indices,2)
    !do j = 1, size(indices,1)
    !     print *, indices(j,1), '->', indices(j,2), weights(j)
    !enddo

    !----------------------------------------------------
    ! Call into online XGrid generation
    !----------------------------------------------------
    xgrid = ESMF_XGridCreate(sideAGrid=(/grid_atm/), sideBGrid=(/grid_ocn/), &
        storeOverlay=.false.,&
        rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    f_tmp = ESMF_FieldCreate(xgrid, xgridside=ESMF_XGRIDSIDE_A, typekind=ESMF_TYPEKIND_R8, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldDestroy(f_tmp)

    ! make sure serialize and deserialize works
    call checkProxy(xgrid, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    ! useful statistics
    call ESMF_XGridGet(xgrid, localDE=0, elementCount=eleCount, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_VMAllReduce(vm, (/eleCount/), totCount, 1, ESMF_REDUCE_SUM, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    if(lpet == 0) write(*, '(A, I4)') 'Total # of XGrid cells: ', totCount(1)

    call ESMF_XGridGet(xgrid, distgridM=distgridM, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    !call ESMF_DistGridPrint(distgridM, rc=localrc)
    !if (ESMF_LogFoundError(localrc, &
    !    ESMF_ERR_PASSTHRU, &
    !    ESMF_CONTEXT, rcToReturn=rc)) return

    ! create a Field on the xgrid
    f_xgrid = ESMF_FieldCreate(xgrid=xgrid, TYPEKIND=ESMF_TYPEKIND_R8, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    !----------------------------------------------------
    ! regrid through the xgrid
    ! set up src flux
    !----------------------------------------------------
    call ESMF_FieldGet(f_atm, farrayPtr=atm, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldGet(f_xgrid, farrayPtr=exf, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldGet(f_ocn, farrayPtr=ocn, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldGet(srcFrac, farrayPtr=srcFracPtr, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldGet(dstFrac, farrayPtr=dstFracPtr, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    !print *, 'src: ', srcFracPtr(:,1)
    !print *, 'dst: ', dstFracPtr

    call ESMF_FieldRegridGetArea(srcArea, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldGet(srcArea, farrayPtr=srcAreaPtr, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_FieldRegridGetArea(dstArea, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldGet(dstArea, farrayPtr=dstAreaPtr, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    mesh_atm = ESMF_GridToMesh(grid_atm, &
      ESMF_STAGGERLOC_CORNER, 0, .false., &
      regridConserve=ESMF_REGRID_CONSERVE_ON, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    mesh_ocn = ESMF_GridToMesh(grid_ocn, &
      ESMF_STAGGERLOC_CORNER, 0, .false., &
      regridConserve=ESMF_REGRID_CONSERVE_ON, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    !----------------------------------------------------
    ! Compute flux integrals
    !----------------------------------------------------
    call init_src(grid_atm, coordX, coordY, atm_nx, atm_ny, atm_dx, atm_dy, atm, rc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    if(lpet == 0) write(*, '(A, E17.10)') 'Total src flux: ', sum(atm)*atm_dx*atm_dy
    !if(atm_nx == 80 .or. (atm_nx == 4 .and. atm_ny == 6)) then
    !  call display_flux2D(atm, srcAreaPtr, srcFracPtr)
    !  call ESMF_FieldWrite(f_atm, file='src.nc', rc=localrc)
    !  if (ESMF_LogFoundError(localrc, &
    !      ESMF_ERR_PASSTHRU, &
    !      ESMF_CONTEXT, rcToReturn=rc)) return
    !  if(lpet == 0) write(*, '(A, E17.10)') 'Total src flux: ', sum(atm)*atm_dx*atm_dy
    !endif                getline(cin, cmd);
    !  error = 0.
    !  write(*, '(A,4I5)') 'bounds: ', lbound(coordX,1), ubound(coordX,1), &
    !    lbound(coordY,1), ubound(coordY,1)
    !  do i = lbound(coordX,1), ubound(coordX,1)
    !    do j = lbound(coordY,1), ubound(coordY,1)
    !      error = error + at                getline(cin, cmd);
    !    enddo
    !  enddo
    !  if(lpet == 0) print *, ' src flux: ', error*atm_dx*atm_dy

    ! X center 
    call ESMF_GridGetCoord(grid_atm, localDE=0, &
        staggerLoc=ESMF_STAGGERLOC_CENTER, &
        coordDim=1, farrayPtr=coordX, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    ! Y center 
    call ESMF_GridGetCoord(grid_atm, localDE=0, &
        staggerLoc=ESMF_STAGGERLOC_CENTER, &
        coordDim=2, farrayPtr=coordY, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    srcsum = 0.
    do i = lbound(coordX,1), ubound(coordX,1)
      do j = lbound(coordY,1), ubound(coordY,1)
        srcsum(1) = srcsum(1) + atm(i,j)*srcAreaPtr(i,j)*srcFracPtr(i,j)
        srcsum(2) = srcsum(2) +          srcAreaPtr(i,j)*srcFracPtr(i,j)
        srcsum(3) = srcsum(3) +          srcAreaPtr(i,j)
      enddo
    enddo
    call ESMF_VMAllReduce(vm, srcsum, allsrcsum, 3, ESMF_REDUCE_SUM, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    if(lpet == 0) print *, ' src flux and area: ', allsrcsum
    call compute_flux2D(vm, atm, srcAreaPtr, srcFracPtr, allsrcsum, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    if(lpet == 0) print *, ' src flux and area sub: ', allsrcsum

    !----------------------------------------------------
    ! call direct regrid
    !----------------------------------------------------
    call ESMF_FieldRegrid(srcField=f_atm, dstField=f_ocn, routehandle=rh_a2o, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call compute_flux2D(vm, ocn, dstAreaPtr, dstFracPtr, allsrcsum, dstflux=.true.,rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    if(lpet == 0) print *, ' dst flux and area from direct Regrid: ', allsrcsum
    dstFlux_reg = allsrcsum(1)
    call ESMF_FieldRegrid(srcField=f_ocn, dstField=f_atm, routehandle=rh_o2a, &
      zeroregion=ESMF_REGION_SELECT, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call compute_flux2D(vm, atm, srcAreaPtr, srcFracPtr, allsrcsum, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    if(lpet == 0) print *, ' recomputed src flux and area from direct Regrid: ', allsrcsum

    !----------------------------------------------------
    ! compute regrid routehandle
    !----------------------------------------------------
    !call ESMF_LogWrite(tag, ESMF_LOGMSG_INFO)
    ! Test passing in srcTermProcessing and pipeLineDepth
    srcTermProcessing=-1 ! Set to -1 to use autotuning
    pipeLineDepth=-1     ! Set to -1 to use autotuning
    call ESMF_FieldRegridStore(xgrid, f_atm, f_xgrid, &
         srcTermProcessing=srcTermProcessing, &
         pipeLineDepth=pipeLineDepth, &
         routehandle=rh_a2x, rc=localrc)
    ! Debug output
    !write(*,*) "srcTermProcessing=",srcTermProcessing
    !write(*,*) "pipeLineDepth=",pipeLineDepth
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldRegridStore(xgrid, f_xgrid, f_atm, routehandle=rh_x2a, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldRegridStore(xgrid, f_ocn, f_xgrid, routehandle=rh_o2x, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldRegridStore(xgrid, f_xgrid, f_ocn, routehandle=rh_x2o, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    !----------------------------------------------------
    ! execute regrid
    ! reinitialize atm data because we called regridstore/regrid on f_atm
    !----------------------------------------------------
    call init_src(grid_atm, coordX, coordY, atm_nx, atm_ny, atm_dx, atm_dy, atm, rc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldRegrid(f_atm, f_xgrid, routehandle=rh_a2x, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    !print *, 'a2x exf', lpet, exf

    allocate(xArea(lbound(exf,1):ubound(exf,1)))
    allocate(xFrac(lbound(exf,1):ubound(exf,1)))
    call ESMF_XGridGet(xgrid, area=xArea, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    xFrac = 1.0

    call compute_flux1D(vm, exf, xArea, xFrac, allsrcsum, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    if(lpet == 0) print *, ' xgrid flux and area: ', allsrcsum
    !call display_flux1D(exf, xArea, xFrac)
    call ESMF_XGridGet(xgrid, sparseMatA2X=sparseMat, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    !do j = 1, size(sparseMat(1)%factorList,1)
    !     print *, sparseMat(1)%factorIndexList(1,j), '->', &
    !      sparseMat(1)%factorIndexList(2,j), sparseMat(1)%factorList(j)
    !enddo
    

    call ESMF_FieldRegrid(f_xgrid, f_ocn, routehandle=rh_x2o, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    !print *, 'x2o ocn', lpet, ocn

    ! Compute flux integrals
    call compute_flux2D(vm, ocn, dstAreaPtr, dstFracPtr, allsrcsum, dstflux=.true.,rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    if(lpet == 0) print *, ' dst flux and area: ', allsrcsum
    dstFlux = allsrcsum(1)

    error = abs(dstFlux - dstFlux_reg)/abs(dstFlux_reg) 
    if(lpet == 0) write(*,'(A,3E18.10)') 'Verify: ', dstFlux_reg, dstFlux, error
    if(error > 1.e-10) then
      print *, 'Regrid through XGrid doesnot agree with direct Regrid'
      if(present(rc)) rc = ESMF_RC_NOT_VALID
      return
    endif

    !call ESMF_FieldWrite(f_atm, 'atm', rc=localrc)
    !if (ESMF_LogFoundError(localrc, &
    !    ESMF_ERR_PASSTHRU, &
    !    ESMF_CONTEXT, rcToReturn=rc)) return
    !call ESMF_FieldWrite(f_ocn, 'ocn', rc=localrc)
    !if (ESMF_LogFoundError(localrc, &
    !    ESMF_ERR_PASSTHRU, &
    !    ESMF_CONTEXT, rcToReturn=rc)) return

    !----------------------------------------------------
    ! Now reverse the regridding direction through xgrid
    !----------------------------------------------------
    !call display_flux1D(exf, xArea, xFrac)
    call ESMF_FieldRegrid(f_ocn, f_xgrid, routehandle=rh_o2x, &
      zeroregion=ESMF_REGION_SELECT, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    !print *, 'o2x exf', lpet, exf
    call compute_flux1D(vm, exf, xArea, xFrac, allsrcsum, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    if(lpet == 0) print *, ' recomputed xgrid flux and area: ', allsrcsum
    !call display_flux1D(exf, xArea, xFrac)

    !call display_flux2D(atm, srcAreaPtr, srcFracPtr)
    call ESMF_FieldRegrid(f_xgrid, f_atm, routehandle=rh_x2a, &
      zeroregion=ESMF_REGION_SELECT, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    !print *, 'x2a atm', lpet, atm
    call compute_flux2D(vm, atm, srcAreaPtr, srcFracPtr, allsrcsum, dstflux=.true.,rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    if(lpet == 0) print *, ' recomputed src flux and area: ', allsrcsum
    !call display_flux2D(atm, srcAreaPtr, srcFracPtr)

    !----------------------------------------------------
    ! clean up
    !----------------------------------------------------
    call ESMF_FieldDestroy(f_atm, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_FieldDestroy(f_ocn, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_GridDestroy(grid_atm, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_GridDestroy(grid_ocn, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_FieldDestroy(f_xgrid, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_XGridDestroy(xgrid, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_MeshDestroy(mesh_atm, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_MeshDestroy(mesh_ocn, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_FieldDestroy(srcFrac, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldDestroy(dstFrac, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

!BOB    deallocate(xArea, xFrac)

    !----------------------------------------------------
    ! release routehandles
    !----------------------------------------------------
    call ESMF_FieldRegridRelease(rh_a2x, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldRegridRelease(rh_x2o, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldRegridRelease(rh_o2x, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldRegridRelease(rh_x2a, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldRegridRelease(rh_a2o, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldRegridRelease(rh_o2a, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    if(present(rc)) rc = ESMF_SUCCESS

  end subroutine test_regrid2xg_online

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

  subroutine test_regrid2xgSph(atm_nx, atm_ny, ocn_nx, ocn_ny, atm_dx, atm_dy, &
    ocn_dx, ocn_dy, atm_sx, atm_sy, ocn_sx, ocn_sy, tag, maxnpet, minnpet, scheme, rc)
    ! arguments
    integer, intent(in)                       :: atm_nx, atm_ny, ocn_nx, ocn_ny
    real(ESMF_KIND_R4), intent(in)            :: atm_dx, atm_dy, ocn_dx, ocn_dy
    real(ESMF_KIND_R4), intent(in), optional  :: atm_sx, atm_sy, ocn_sx, ocn_sy
    character(len=*), intent(in), optional    :: tag
    integer, intent(in) , optional            :: maxnpet, minnpet
    integer, intent(in) , optional            :: scheme
    integer, intent(out), optional            :: rc

    ! local variables
    type(ESMF_Grid)                 :: grid_atm, grid_ocn
    type(ESMF_Field)                :: f_atm, f_ocn, f_xgrid, f_error, f_ocn_dir
    real(ESMF_KIND_R8)              :: startx, starty
    integer                         :: localrc, npet, i, j, lpet
    real(ESMF_KIND_R8), pointer     :: weights(:)
    integer(ESMF_KIND_I4), pointer  :: indices(:,:)
    real(ESMF_KIND_R8), pointer     :: coordX(:,:), coordY(:,:)
    type(ESMF_XGrid)                :: xgrid

    type(ESMF_XGridSpec)            :: sparseMatA2X(1)
    integer                         :: gn(2), simax, dimax, l_minnpet, l_maxnpet
    real(ESMF_KIND_R8), pointer     :: atm(:,:), ocn(:,:), ptr_error(:,:), exf(:), xArea(:), xFrac(:)
    real(ESMF_KIND_R8), pointer     :: ocn_dir(:,:)
    real(ESMF_KIND_R8), pointer     :: srcFracPtr(:,:), dstFracPtr(:,:)
    real(ESMF_KIND_R8), pointer     :: srcAreaPtr(:,:), dstAreaPtr(:,:)
    type(ESMF_RouteHandle)          :: rh_a2o, rh_o2a, rh_a2x, rh_x2o, rh_o2x, rh_x2a
    type(ESMF_DistGrid)             :: distgridM, distgrid1, distgrid2
    type(ESMF_XGridSpec)            :: sparseMat(1)

    type(ESMF_Field)                :: fa_atm, fa_ocn, fa_xgrid
    type(ESMF_Field)                :: aa_atm, aa_ocn, aa_xgrid, fieldM
    type(ESMF_Field)                :: srcFrac, dstFrac, srcArea, dstArea
    type(ESMF_Mesh)                 :: mesh_atm, mesh_ocn, mesh_xgrid, meshM
    real(ESMF_KIND_R8)              :: srcsum(3), allsrcsum(3), scale
    real(ESMF_KIND_R8)              :: dstFlux_reg, dstFlux, totalXArea, totalSrcArea, error
    integer                         :: l_scheme, eleCount

    type(ESMF_VM)   :: vm

    localrc = ESMF_SUCCESS
    if(present(rc)) rc = ESMF_SUCCESS

    call ESMF_VMGetCurrent(vm, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_VMGet(vm, petCount=npet, localPet=lpet, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    ! Don't run this test if there isn't enough pets or too many cuts
    l_minnpet = 1
    l_maxnpet = 128
    if(present(minnpet)) l_minnpet = minnpet
    if(present(maxnpet)) l_maxnpet = maxnpet
    if(npet > l_maxnpet .or. npet < l_minnpet) return

    if(lpet == 0 .and. present(tag)) then
      print *, '!------------------------------------'
      print *, '!', tag
      print *, '!------------------------------------'
    endif

    l_scheme = ESMF_REGRID_SCHEME_REGION3D
    if(present(scheme)) l_scheme = scheme

    ! Grid constants
    ! Atmosphere covers the area (-165,30) - (-15, 50)
    ! North Pacific Ocean covers (-165,30) - (-120, 50)
    ! Running the Grids from NW corner to SE corner, 
    ! change the starting Y coord and dy
    !
    ! 1 degree atmosphere, 1 degree ocean
    !atm_nx = 5
    !atm_ny = 20
    !atm_dx = 2.
    !atm_dy = 1.
    !ocn_nx = 5
    !ocn_ny = 20
    !ocn_dx = 1.
    !ocn_dy = 1.
    !
    !atm_sx = -165.
    !atm_sy = 30.
    !ocn_sx = -165.
    !ocn_sy = 30.
    
    !------------- ATM --------------
    ! atm grid, horizontally decomposed
    if(l_scheme == ESMF_REGRID_SCHEME_FULL3D) then
      grid_atm = ESMF_GridCreate1PeriDim(maxIndex=(/atm_nx, atm_ny/), &
        indexflag=ESMF_INDEX_GLOBAL, &
        gridEdgeLWidth=(/0,0/), gridEdgeUWidth=(/0,1/), &
        regDecomp=(/npet, 1/), &
        rc=localrc)
    else
      grid_atm = ESMF_GridCreateNoPeriDim(maxIndex=(/atm_nx, atm_ny/), &
        indexflag=ESMF_INDEX_GLOBAL, &
        gridEdgeLWidth=(/0,0/), gridEdgeUWidth=(/1,1/), &
        regDecomp=(/npet, 1/), &
        rc=localrc)
    endif 
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_GridAddCoord(grid_atm, staggerloc=ESMF_STAGGERLOC_CENTER, &
        rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_GridAddCoord(grid_atm, staggerloc=ESMF_STAGGERLOC_CORNER, &
        rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    ! global indexing
    ! atm grid is not decomposed in the y direction
    !startx = lpet*atm_nx/npet*atm_dx
    startx = atm_sx
    starty = atm_sy
    ! compute coord
    ! X center
    call ESMF_GridGetCoord(grid_atm, localDE=0, staggerLoc=ESMF_STAGGERLOC_CENTER, &
        coordDim=1, farrayPtr=coordX, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    ! Y center
    call ESMF_GridGetCoord(grid_atm, localDE=0, staggerLoc=ESMF_STAGGERLOC_CENTER, &
        coordDim=2, farrayPtr=coordY, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    do i = lbound(coordX,1), ubound(coordX,1)
      do j = lbound(coordX, 2), ubound(coordX, 2)
        coordX(i,j) = startx + atm_dx/2. + (i-1)*atm_dx
        coordY(i,j) = starty + atm_dy/2. + (j-1)*atm_dy
      enddo
    enddo
    print *, 'startx: ', startx, lbound(coordX, 1), ubound(coordX, 1), 'coordX: ', coordX(:,1)
    ! X corner
    call ESMF_GridGetCoord(grid_atm, localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, &
        coordDim=1, farrayPtr=coordX, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    ! Y corner
    call ESMF_GridGetCoord(grid_atm, localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, &
        coordDim=2, farrayPtr=coordY, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    do i = lbound(coordX,1), ubound(coordX,1)
      do j = lbound(coordX, 2), ubound(coordX, 2)
        coordX(i,j) = startx + (i-1)*atm_dx
        coordY(i,j) = starty + (j-1)*atm_dy
      enddo
    enddo

    !------------- OCN --------------
    ! ocn grid, horizontally decomposed
    if(l_scheme == ESMF_REGRID_SCHEME_FULL3D) then
      grid_ocn = ESMF_GridCreate1PeriDim(maxIndex=(/ocn_nx, ocn_ny/), &
        indexflag=ESMF_INDEX_GLOBAL, &
        gridEdgeLWidth=(/0,0/), gridEdgeUWidth=(/0,1/), &
        regDecomp=(/npet, 1/), &
        rc=localrc)
    else
      grid_ocn = ESMF_GridCreateNoPeriDim(maxIndex=(/ocn_nx, ocn_ny/), &
        indexflag=ESMF_INDEX_GLOBAL, &
        gridEdgeLWidth=(/0,0/), gridEdgeUWidth=(/1,1/), &
        regDecomp=(/npet, 1/), &
        rc=localrc)
    endif 
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_GridAddCoord(grid_ocn, staggerloc=ESMF_STAGGERLOC_CENTER, &
        rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_GridAddCoord(grid_ocn, staggerloc=ESMF_STAGGERLOC_CORNER, &
        rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    ! ocn grid is not decomposed in the y direction
    !startx = lpet*ocn_nx/npet*ocn_dx
    startx = ocn_sx
    starty = ocn_sy
    ! compute coord
    ! X center
    call ESMF_GridGetCoord(grid_ocn, localDE=0, staggerLoc=ESMF_STAGGERLOC_CENTER, &
        coordDim=1, farrayPtr=coordX, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    ! Y center
    call ESMF_GridGetCoord(grid_ocn, localDE=0, staggerLoc=ESMF_STAGGERLOC_CENTER, &
        coordDim=2, farrayPtr=coordY, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    do i = lbound(coordX,1), ubound(coordX,1)
      do j = lbound(coordX, 2), ubound(coordX, 2)
        coordX(i,j) = startx + ocn_dx/2. + (i-1)*ocn_dx
        coordY(i,j) = starty + ocn_dy/2. + (j-1)*ocn_dy
      enddo
    enddo
    !print *, 'startx: ', startx, lbound(coordX, 1), 'coordX: ', coordX(:,1), coordX(1,:)
    ! X corner
    call ESMF_GridGetCoord(grid_ocn, localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, &
        coordDim=1, farrayPtr=coordX, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    ! Y corner
    call ESMF_GridGetCoord(grid_ocn, localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, &
        coordDim=2, farrayPtr=coordY, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    do i = lbound(coordX,1), ubound(coordX,1)
      do j = lbound(coordX, 2), ubound(coordX, 2)
        coordX(i,j) = startx + (i-1)*ocn_dx
        coordY(i,j) = starty + (j-1)*ocn_dy
      enddo
    enddo

    ! build Fields on the Grids
    f_atm = ESMF_FieldCreate(grid_atm, typekind=ESMF_TYPEKIND_R8, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    f_error = ESMF_FieldCreate(grid_atm, typekind=ESMF_TYPEKIND_R8, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    f_ocn = ESMF_FieldCreate(grid_ocn, typekind=ESMF_TYPEKIND_R8, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    f_ocn_dir = ESMF_FieldCreate(grid_ocn, typekind=ESMF_TYPEKIND_R8, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    srcFrac = ESMF_FieldCreate(grid_atm, typekind=ESMF_TYPEKIND_R8, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    dstFrac = ESMF_FieldCreate(grid_ocn, typekind=ESMF_TYPEKIND_R8, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    srcArea = ESMF_FieldCreate(grid_atm, typekind=ESMF_TYPEKIND_R8, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    dstArea = ESMF_FieldCreate(grid_ocn, typekind=ESMF_TYPEKIND_R8, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    !----------------------------------------------------
    ! Call RegridStore here to compute SMM weights and indices and direct RHs
    !----------------------------------------------------

    call ESMF_FieldRegridStore(srcField=f_atm, dstField=f_ocn, &
      routehandle=rh_a2o, &
      regridmethod=ESMF_REGRIDMETHOD_CONSERVE, &
      unmappedaction = ESMF_UNMAPPEDACTION_IGNORE, &
      srcFracField=srcFrac, dstFracField=dstFrac, & 
      factorIndexList=indices, factorList=weights, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    ! make sure the numbers are consistent
    !print *, lpet, 'weights', size(weights)
    !print *, lpet, 'indices', size(indices,1),'-', size(indices,2)
    !do j = 1, size(indices,1)
    !     print *, indices(j,1), '->', indices(j,2), weights(j)
    !enddo

    call ESMF_FieldRegridStore(srcField=f_ocn, dstField=f_atm, &
      routehandle=rh_o2a, &
      regridmethod=ESMF_REGRIDMETHOD_CONSERVE, &
      unmappedaction = ESMF_UNMAPPEDACTION_IGNORE, &
      factorIndexList=indices, factorList=weights, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    ! make sure the numbers are consistent
    !print *, lpet, 'weights', size(weights)
    !print *, lpet, 'indices', size(indices,1),'-', size(indices,2)
    !do j = 1, size(indices,1)
    !     print *, indices(j,1), '->', indices(j,2), weights(j)
    !enddo

    !----------------------------------------------------
    ! Call into online generation
    !----------------------------------------------------
    xgrid = ESMF_XGridCreate(sideAGrid=(/grid_atm/), sideBGrid=(/grid_ocn/), &
        storeOverlay=.true.,&
        rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call checkProxy(xgrid, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    ! Query the XGrid for the overlay mesh
    call ESMF_XGridGet(xgrid, mesh=meshM, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    ! Create a field on the overlay mesh
    fieldM = ESMF_FieldCreate(mesh=meshM, typekind=ESMF_TYPEKIND_R8, &
        rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_XGridGet(xgrid, distgridM=distgridM, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    !call ESMF_DistGridPrint(distgridM, rc=localrc)
    !if (ESMF_LogFoundError(localrc, &
    !    ESMF_ERR_PASSTHRU, &
    !    ESMF_CONTEXT, rcToReturn=rc)) return

    ! create a Field on the xgrid
    f_xgrid = ESMF_FieldCreate(xgrid=xgrid, TYPEKIND=ESMF_TYPEKIND_R8, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    !----------------------------------------------------
    ! regrid through the xgrid
    !----------------------------------------------------

    !----------------------------------------------------
    ! set up src flux
    !----------------------------------------------------
    call ESMF_FieldGet(f_atm, farrayPtr=atm, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldGet(f_error, farrayPtr=ptr_error, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldGet(f_xgrid, farrayPtr=exf, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldGet(f_ocn, farrayPtr=ocn, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldGet(f_ocn_dir, farrayPtr=ocn_dir, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldGet(srcFrac, farrayPtr=srcFracPtr, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldGet(dstFrac, farrayPtr=dstFracPtr, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    !print *, 'src: ', srcFracPtr(:,1)
    !print *, 'dst: ', dstFracPtr

    call ESMF_XGridGet(xgrid, elementCount=eleCount, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    write(*, *) 'eleCount = ', eleCount, 'size(fptr) = ', size(exf)
    if(eleCount /= size(exf)) then
      print *, 'elementCount is not equal to size(exf)'
      if(present(rc)) rc = ESMF_RC_NOT_VALID
      return
    endif

    call ESMF_FieldRegridGetArea(srcArea, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldGet(srcArea, farrayPtr=srcAreaPtr, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_FieldRegridGetArea(dstArea, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldGet(dstArea, farrayPtr=dstAreaPtr, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    mesh_atm = ESMF_GridToMesh(grid_atm, &
      ESMF_STAGGERLOC_CORNER, 0, .false., &
      regridConserve=ESMF_REGRID_CONSERVE_ON, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    mesh_ocn = ESMF_GridToMesh(grid_ocn, &
      ESMF_STAGGERLOC_CORNER, 0, .false., &
      regridConserve=ESMF_REGRID_CONSERVE_ON, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    !----------------------------------------------------
    ! Compute flux integrals
    !----------------------------------------------------
    call init_src_sph(grid_atm, coordX, coordY, atm_nx, atm_ny, atm_dx, atm_dy, atm, rc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_GridGetCoord(grid_atm, localDE=0, &
        staggerLoc=ESMF_STAGGERLOC_CENTER, &
        coordDim=1, farrayPtr=coordX, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    ! Y center 
    call ESMF_GridGetCoord(grid_atm, localDE=0, &
        staggerLoc=ESMF_STAGGERLOC_CENTER, &
        coordDim=2, farrayPtr=coordY, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    srcsum = 0.
    do i = lbound(coordX,1), ubound(coordX,1)
      do j = lbound(coordX,2), ubound(coordX,2)
        !atm(i,j) = 0.5*sin(coordX(i)/(atm_nx*atm_dx)*2*3.14159)*cos(coordY(j)&
        !/(atm_ny*atm_dy)*4*3.14159)
        !atm(i,j) = i
        srcsum(1) = srcsum(1) + atm(i,j)*srcAreaPtr(i,j)*srcFracPtr(i,j)
        srcsum(2) = srcsum(2) +          srcAreaPtr(i,j)*srcFracPtr(i,j)
        srcsum(3) = srcsum(3) +          srcAreaPtr(i,j)
      enddo
    enddo
    call ESMF_VMAllReduce(vm, srcsum, allsrcsum, 3, ESMF_REDUCE_SUM, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    !print *, lpet, ' src flux and area: ', atm
    scale = 3.1415926535897932/180.
    totalSrcArea = allsrcsum(3)
    if(lpet == 0) print *, ' src flux and area: ', allsrcsum, &
      (sin(scale*(atm_sy+atm_ny*atm_dy))-sin(scale*atm_sy))*(atm_nx*atm_dx)*scale
    call compute_flux2D(vm, atm, srcAreaPtr, srcFracPtr, allsrcsum, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    if(lpet == 0) print *, ' src flux and area sub: ', allsrcsum

    !----------------------------------------------------
    ! call direct regrid
    !----------------------------------------------------
    !call display_flux2D(atm, srcAreaPtr, srcFracPtr)
    call ESMF_FieldRegrid(srcField=f_atm, dstField=f_ocn, routehandle=rh_a2o, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call compute_flux2D(vm, ocn, dstAreaPtr, dstFracPtr, allsrcsum, dstflux=.true., rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    if(lpet == 0) print *, ' dst flux and area from direct Regrid: ', allsrcsum
    dstFlux_reg = allsrcsum(1)
    ocn_dir = ocn

    !call display_flux2D(ocn, dstAreaPtr, dstFracPtr)
    call ESMF_FieldRegrid(srcField=f_ocn, dstField=f_atm, routehandle=rh_o2a, &
      zeroregion=ESMF_REGION_SELECT, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call compute_flux2D(vm, atm, srcAreaPtr, srcFracPtr, allsrcsum, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    if(lpet == 0) print *, ' recomputed src flux and area from direct Regrid: ', allsrcsum

    call init_src_sph(grid_atm, coordX, coordY, atm_nx, atm_ny, atm_dx, atm_dy, ptr_error, rc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call diff_ptr_sph(grid_atm, coordX, coordY, atm_nx, atm_ny, atm_dx, atm_dy, atm, ptr_error, rc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    !if(l_scheme == ESMF_REGRID_SCHEME_FULL3D) then
    !  call ESMF_FieldWrite(f_error, file='error.nc', rc=localrc)
    !  if (ESMF_LogFoundError(localrc, &
    !      ESMF_ERR_PASSTHRU, &
    !      ESMF_CONTEXT, rcToReturn=rc)) return
    !endif

    !----------------------------------------------------
    ! compute xgrid regrid routehandle
    !----------------------------------------------------
    call ESMF_FieldRegridStore(xgrid, f_atm, f_xgrid, routehandle=rh_a2x, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldRegridStore(xgrid, f_xgrid, f_atm, routehandle=rh_x2a, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldRegridStore(xgrid, f_ocn, f_xgrid, routehandle=rh_o2x, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldRegridStore(xgrid, f_xgrid, f_ocn, routehandle=rh_x2o, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    !----------------------------------------------------
    ! execute regrid
    ! reinitialize atm data because we called regridstore/regrid on f_atm
    !----------------------------------------------------
    call init_src_sph(grid_atm, coordX, coordY, atm_nx, atm_ny, atm_dx, atm_dy, atm, rc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldRegrid(f_atm, f_xgrid, routehandle=rh_a2x, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    !print *, 'a2x exf', lpet, exf

    allocate(xArea(lbound(exf,1):ubound(exf,1)))
    allocate(xFrac(lbound(exf,1):ubound(exf,1)))
    call ESMF_XGridGet(xgrid, area=xArea, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    xFrac = 1.0
    call compute_flux1D(vm, exf, xArea, xFrac, allsrcsum, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    totalXArea = allsrcsum(3)
    if(lpet == 0) print *, ' xgrid flux and area: ', allsrcsum
    !call display_flux1D(exf, xArea, xFrac)
    call ESMF_XGridGet(xgrid, sparseMatA2X=sparseMat, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    !do j = 1, size(sparseMat(1)%factorList,1)
    !     print *, sparseMat(1)%factorIndexList(1,j), '->', &
    !      sparseMat(1)%factorIndexList(2,j), sparseMat(1)%factorList(j)
    !enddo
    

    call ESMF_FieldRegrid(f_xgrid, f_ocn, routehandle=rh_x2o, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    !print *, 'x2o ocn', lpet, ocn

    ! Compute flux integrals
    call compute_flux2D(vm, ocn, dstAreaPtr, dstFracPtr, allsrcsum, dstflux=.true., rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    if(lpet == 0) print *, ' dst flux and area: ', allsrcsum
    dstFlux = allsrcsum(1)

    !----------------------------------------------------
    ! Verify
    !----------------------------------------------------
    call diff_ptr_sph(grid_ocn, coordX, coordY, ocn_nx, ocn_ny, ocn_dx, ocn_dy, ocn, ocn_dir, rc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    !if(l_scheme == ESMF_REGRID_SCHEME_FULL3D) then
    !  call ESMF_FieldWrite(f_ocn_dir, file='xgrid_regrid_diff.nc', rc=localrc)
    !  if (ESMF_LogFoundError(localrc, &
    !      ESMF_ERR_PASSTHRU, &
    !      ESMF_CONTEXT, rcToReturn=rc)) return
    !endif
    error = abs(dstFlux - dstFlux_reg)/abs(dstFlux_reg) 
    if(lpet == 0) write(*,'(A,4E18.10)') 'Verify: ', dstFlux_reg, dstFlux, error, abs(totalXArea-totalSrcArea)/totalSrcArea
    if(error > 1.e-4) then
      print *, 'Regrid through XGrid doesnot agree with direct Regrid'
      if(present(rc)) rc = ESMF_RC_NOT_VALID
      return
    endif

    !call ESMF_FieldWrite(f_atm, 'atm', rc=localrc)
    !if (ESMF_LogFoundError(localrc, &
    !    ESMF_ERR_PASSTHRU, &
    !    ESMF_CONTEXT, rcToReturn=rc)) return
    !call ESMF_FieldWrite(f_ocn, 'ocn', rc=localrc)
    !if (ESMF_LogFoundError(localrc, &
    !    ESMF_ERR_PASSTHRU, &
    !    ESMF_CONTEXT, rcToReturn=rc)) return

    !----------------------------------------------------
    ! Now reverse the regridding direction through xgrid
    !----------------------------------------------------
    !call display_flux1D(exf, xArea, xFrac)
    call ESMF_FieldRegrid(f_ocn, f_xgrid, routehandle=rh_o2x, &
      zeroregion=ESMF_REGION_SELECT, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    !print *, 'o2x exf', lpet, exf
    call compute_flux1D(vm, exf, xArea, xFrac, allsrcsum, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    if(lpet == 0) print *, ' recomputed xgrid flux and area: ', allsrcsum
    !call display_flux1D(exf, xArea, xFrac)

    !call display_flux2D(atm, srcAreaPtr, srcFracPtr)
    call ESMF_FieldRegrid(f_xgrid, f_atm, routehandle=rh_x2a, &
      zeroregion=ESMF_REGION_SELECT, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    !print *, 'x2a atm', lpet, atm
    call compute_flux2D(vm, atm, srcAreaPtr, srcFracPtr, allsrcsum, dstflux=.true., rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    if(lpet == 0) print *, ' recomputed src flux and area: ', allsrcsum
    !call display_flux2D(atm, srcAreaPtr, srcFracPtr)

    call init_src_sph(grid_atm, coordX, coordY, atm_nx, atm_ny, atm_dx, atm_dy, ptr_error, rc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call diff_ptr_sph(grid_atm, coordX, coordY, atm_nx, atm_ny, atm_dx, atm_dy, atm, ptr_error, rc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    !if(l_scheme == ESMF_REGRID_SCHEME_FULL3D) then
    !  call ESMF_FieldWrite(f_error, file='error1.nc', rc=localrc)
    !  if (ESMF_LogFoundError(localrc, &
    !      ESMF_ERR_PASSTHRU, &
    !      ESMF_CONTEXT, rcToReturn=rc)) return
    !endif

    !----------------------------------------------------
    ! clean up
    !----------------------------------------------------
    call ESMF_FieldDestroy(f_atm, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_FieldDestroy(f_ocn, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_GridDestroy(grid_atm, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_GridDestroy(grid_ocn, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_FieldDestroy(f_xgrid, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_XGridDestroy(xgrid, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_MeshDestroy(mesh_atm, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_MeshDestroy(mesh_ocn, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_FieldDestroy(srcFrac, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldDestroy(dstFrac, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

!BOB    deallocate(xArea, xFrac)

    !----------------------------------------------------
    ! release routehandles
    !----------------------------------------------------
    call ESMF_FieldRegridRelease(rh_a2x, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldRegridRelease(rh_x2o, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldRegridRelease(rh_o2x, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldRegridRelease(rh_x2a, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldRegridRelease(rh_a2o, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    call ESMF_FieldRegridRelease(rh_o2a, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    if(present(rc)) rc = ESMF_SUCCESS

  end subroutine test_regrid2xgSph

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

  subroutine compute_flux1D(vm, flux_density, area, fraction, allsum, rc)
    type(ESMF_VM), intent(in)        :: vm
    real(ESMF_KIND_R8), pointer      :: flux_density(:) 
    real(ESMF_KIND_R8), pointer      :: area(:) 
    real(ESMF_KIND_R8), pointer      :: fraction(:) 
    real(ESMF_KIND_R8), intent(out)  :: allsum(3)
    integer, intent(out), optional   :: rc

    real(ESMF_KIND_R8)               :: sum(3)
    integer                          :: i,j, localrc

    if(present(rc)) rc = ESMF_SUCCESS

    sum = 0.
    do i = lbound(flux_density, 1), ubound(flux_density, 1)
      sum(1) = sum(1) + flux_density(i)*area(i)*fraction(i)
      sum(2) = sum(2) +                 area(i)*fraction(i)
      sum(3) = sum(3) +                 area(i)
    enddo

    call ESMF_VMAllReduce(vm, sum, allsum, 3, ESMF_REDUCE_SUM, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

  end subroutine compute_flux1D

  subroutine compute_flux2D(vm, flux_density, area, fraction, allsum, dstflux, rc)
    type(ESMF_VM), intent(in)        :: vm
    real(ESMF_KIND_R8), pointer      :: flux_density(:,:) 
    real(ESMF_KIND_R8), pointer      :: area(:,:) 
    real(ESMF_KIND_R8), pointer      :: fraction(:,:) 
    real(ESMF_KIND_R8), intent(out)  :: allsum(3)
    logical, intent(in) , optional   :: dstflux
    integer, intent(out), optional   :: rc

    real(ESMF_KIND_R8)               :: sum(3)
    integer                          :: i,j, localrc, npet, lpet
    logical                          :: l_dstflux

    if(present(rc)) rc = ESMF_SUCCESS
    l_dstflux = .false.
    if(present(dstflux)) l_dstflux = dstflux

    call ESMF_VMGet(vm, petCount=npet, localPet=lpet, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
  
    !if(lpet == 0) write(*, '(A, 4I5)') 'compute_flux2D bounds: ', &
    !  lbound(flux_density, 1), ubound(flux_density, 1), &
    !  lbound(flux_density, 2), ubound(flux_density, 2)

    sum = 0.
    do i = lbound(flux_density, 1), ubound(flux_density, 1)
      do j = lbound(flux_density, 2), ubound(flux_density, 2)
        if(l_dstflux) then
          sum(1) = sum(1) + flux_density(i,j)*area(i,j)
        else
          sum(1) = sum(1) + flux_density(i,j)*area(i,j)*fraction(i,j)
        endif
        sum(2) = sum(2) +                 area(i,j)*fraction(i,j)
        sum(3) = sum(3) +                 area(i,j)
      enddo
    enddo

    call ESMF_VMAllReduce(vm, sum, allsum, 3, ESMF_REDUCE_SUM, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

  end subroutine compute_flux2D

  subroutine display_flux1D(flux_density, area, fraction, rc)
    real(ESMF_KIND_R8), pointer      :: flux_density(:) 
    real(ESMF_KIND_R8), pointer      :: area(:) 
    real(ESMF_KIND_R8), pointer      :: fraction(:) 
    integer, intent(out), optional   :: rc

    integer                          :: i,localrc, npet, lpet
    type(ESMF_VM)                    :: vm

    localrc = ESMF_SUCCESS
    if(present(rc)) rc = ESMF_SUCCESS

    call ESMF_VMGetCurrent(vm, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_VMGet(vm, petCount=npet, localPet=lpet, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    do i = lbound(flux_density, 1), ubound(flux_density, 1)
      !print *, i,flux_density(i), area(i) !, fraction(i)
      write(*,'(I3,I7,2F14.10)') lpet,i,flux_density(i), area(i)
    enddo

  end subroutine display_flux1D

  subroutine display_flux2D(flux_density, area, fraction, rc)
    real(ESMF_KIND_R8), pointer      :: flux_density(:,:) 
    real(ESMF_KIND_R8), pointer      :: area(:,:) 
    real(ESMF_KIND_R8), pointer      :: fraction(:,:) 
    integer, intent(out), optional   :: rc

    integer                          :: i,j, localrc, npet, lpet
    type(ESMF_VM)                    :: vm

    localrc = ESMF_SUCCESS
    if(present(rc)) rc = ESMF_SUCCESS

    call ESMF_VMGetCurrent(vm, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_VMGet(vm, petCount=npet, localPet=lpet, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
  
    !if(lpet == 0) write(*, '(A, 4I5)') 'display_flux2D bounds: ', &
    !  lbound(flux_density, 1), ubound(flux_density, 1), &
    !  lbound(flux_density, 2), ubound(flux_density, 2)

    do i = lbound(flux_density, 1), ubound(flux_density, 1)
      do j = lbound(flux_density, 2), ubound(flux_density, 2)
        write(*,'(I3,2I5,3F14.10)') lpet, i,j,flux_density(i,j), area(i,j), fraction(i,j)
      enddo
    enddo

  end subroutine display_flux2D

  subroutine init_src(grid, coordX, coordY, nx, ny, dx, dy, src, rc)
    type(ESMF_Grid), intent(in)      :: grid
    real(ESMF_KIND_R8), pointer      :: coordX(:)
    real(ESMF_KIND_R8), pointer      :: coordY(:)
    integer, intent(in)              :: nx, ny
    real, intent(in)                 :: dx, dy
    real(ESMF_KIND_R8), pointer      :: src(:,:)
    integer, intent(out), optional   :: rc

    integer                          :: localrc, i,j 

    if(present(rc)) rc = ESMF_SUCCESS

    ! X center 
    call ESMF_GridGetCoord(grid, localDE=0, &
        staggerLoc=ESMF_STAGGERLOC_CENTER, &
        coordDim=1, farrayPtr=coordX, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    ! Y center 
    call ESMF_GridGetCoord(grid, localDE=0, &
        staggerLoc=ESMF_STAGGERLOC_CENTER, &
        coordDim=2, farrayPtr=coordY, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    do i = lbound(coordX,1), ubound(coordX,1)
      do j = lbound(coordY,1), ubound(coordY,1)
        src(i,j) = 0.5*sin(coordX(i)/(nx*dx)*3*3.14159)*cos(coordY(j)&
        /(ny*dy)*5*3.14159)+0.5
        !src(i,j) = 2.
      enddo
    enddo

  end subroutine init_src

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

  subroutine init_src_sph(grid, coordX, coordY, nx, ny, dx, dy, src, rc)
    type(ESMF_Grid), intent(in)      :: grid
    real(ESMF_KIND_R8), pointer      :: coordX(:,:)
    real(ESMF_KIND_R8), pointer      :: coordY(:,:)
    integer, intent(in)              :: nx, ny
    real, intent(in)                 :: dx, dy
    real(ESMF_KIND_R8), pointer      :: src(:,:)
    integer, intent(out), optional   :: rc

    integer                          :: localrc, i,j,npet, lpet
    type(ESMF_VM)                    :: vm

    if(present(rc)) rc = ESMF_SUCCESS

    call ESMF_VMGetCurrent(vm, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_VMGet(vm, petCount=npet, localPet=lpet, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    ! X center 
    call ESMF_GridGetCoord(grid, localDE=0, &
        staggerLoc=ESMF_STAGGERLOC_CENTER, &
        coordDim=1, farrayPtr=coordX, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    ! Y center 
    call ESMF_GridGetCoord(grid, localDE=0, &
        staggerLoc=ESMF_STAGGERLOC_CENTER, &
        coordDim=2, farrayPtr=coordY, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    write(*, '(A, 5I5)') 'init_src_sph bounds: ', lpet, &
      lbound(src, 1), ubound(src, 1), &
      lbound(src, 2), ubound(src, 2)

    do i = lbound(coordX,1), ubound(coordX,1)
      do j = lbound(coordX,2), ubound(coordX,2)
        !src(i,j) = 0.5*sin(coordX(i,j)/(nx*dx)*3*3.14159)*cos(coordY(i,j)&
        !/(ny*dy)*5*3.14159)+10.5
        src(i,j) = 2.
      enddo
    enddo

  end subroutine init_src_sph

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

  subroutine diff_ptr_sph(grid, coordX, coordY, nx, ny, dx, dy, src, error, rc)
    type(ESMF_Grid), intent(in)      :: grid
    real(ESMF_KIND_R8), pointer      :: coordX(:,:)
    real(ESMF_KIND_R8), pointer      :: coordY(:,:)
    integer, intent(in)              :: nx, ny
    real, intent(in)                 :: dx, dy
    real(ESMF_KIND_R8), pointer      :: src(:,:), error(:,:)
    integer, intent(out), optional   :: rc

    integer                          :: localrc, i,j 

    if(present(rc)) rc = ESMF_SUCCESS

    call ESMF_GridGetCoord(grid, localDE=0, &
        staggerLoc=ESMF_STAGGERLOC_CENTER, &
        coordDim=1, farrayPtr=coordX, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    ! Y center 
    call ESMF_GridGetCoord(grid, localDE=0, &
        staggerLoc=ESMF_STAGGERLOC_CENTER, &
        coordDim=2, farrayPtr=coordY, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    do i = lbound(coordX,1), ubound(coordX,1)
      do j = lbound(coordX,2), ubound(coordX,2)
        error(i,j) = src(i,j) - error(i,j)
      enddo
    enddo

  end subroutine diff_ptr_sph

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

  subroutine checkProxy(xgrid, rc)
    type(ESMF_XGrid), intent(in)    :: xgrid
    integer, intent(out), optional  :: rc
    
    type(ESMF_XGrid)                :: xgrid1
    character, pointer              :: buffer(:)
    integer                         :: buff_length, offset, localrc

    if(present(rc)) rc = ESMF_SUCCESS
  
    ! Allocate serialization buffer

    buff_length = 1
    allocate (buffer(0:buff_length-1))
    offset = 0
    call ESMF_XGridSerialize(xgrid, buffer, buff_length, offset, &
        inquireflag=ESMF_INQUIREONLY, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return
    deallocate (buffer)

    buff_length = offset
    allocate (buffer(0:buff_length-1))

    ! call serialize and deserialize and verify again

    offset = 0
    call ESMF_XGridSerialize(xgrid, buffer, buff_length, offset, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    offset = 0

    xgrid1 = ESMF_XGridDeserialize(buffer, offset, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    deallocate (buffer)

    call ESMF_XGridValidate(xgrid1, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

  end subroutine checkProxy

end program ESMF_FieldRegridXGUTest