test_regrid2xg_online Subroutine

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

Type IntentOptional Attributes Name
integer, intent(in) :: atm_nx
integer, intent(in) :: atm_ny
integer, intent(in) :: ocn_nx
integer, intent(in) :: ocn_ny
real(kind=ESMF_KIND_R4), intent(in) :: atm_dx
real(kind=ESMF_KIND_R4), intent(in) :: atm_dy
real(kind=ESMF_KIND_R4), intent(in) :: ocn_dx
real(kind=ESMF_KIND_R4), intent(in) :: ocn_dy
real(kind=ESMF_KIND_R4), intent(in), optional :: atm_sx
real(kind=ESMF_KIND_R4), intent(in), optional :: atm_sy
real(kind=ESMF_KIND_R4), intent(in), optional :: ocn_sx
real(kind=ESMF_KIND_R4), intent(in), optional :: ocn_sy
character(len=*), intent(in), optional :: tag
integer, intent(in), optional :: maxnpet
integer, intent(in), optional :: minnpet
type(ESMF_Index_Flag), intent(in), optional :: indexflag
integer, intent(out), optional :: rc

Source Code

  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