ESMF_XGridEx Program

    program ESMF_XGridEx

!------------------------------------------------------------------------------
!ESMF_MULTI_PROC_EXAMPLE        String used by test script to count examples.
!==============================================================================
! !PROGRAM: ESMF_XGridEx - XGrid Examples
!
! !DESCRIPTION:
!
! This program shows examples of XGrid creation and usage
!-----------------------------------------------------------------------------
#include "ESMF.h"
    ! ESMF Framework module
    use ESMF_TestMod
    use ESMF
    use ESMF_XGridMod
    use ESMF_XGridCreateMod
    implicit none

    ! Local variables
    integer                             :: localrc, i, j
    type(ESMF_XGrid)                    :: xgrid
    type(ESMF_Grid)                     :: sideA(2), sideB(1)
    type(ESMF_DistGrid)                 :: sideAdg(2), sideBdg(1), distgrid
    real(ESMF_KIND_R8)                  :: centroid(12,2), area(12)
    type(ESMF_XGridSpec)                :: sparseMatA2X(2), sparseMatX2B(1)

    type(ESMF_Grid)                     :: l_sideA(2), l_sideB(1)
    type(ESMF_DistGrid)                 :: l_sideAdg(2), l_sideBdg(1)
    type(ESMF_XGridSpec)                :: l_sparseMatA2X(2), l_sparseMatX2B(1)
    type(ESMF_Field)                    :: field, srcField(2), dstField(1)
    type(ESMF_Field)                    :: dstFrac, dstFrac2

    integer                             :: eleCount, ngridA, ngridB
    integer                             :: elb, eub, ec

    real(ESMF_KIND_R8), pointer         :: farrayPtr(:,:), xfarrayPtr(:)
    integer                             :: xlb(1), xub(1)
    type(ESMF_RouteHandle)              :: rh_src2xgrid(2), rh_xgrid2dst(1)

    real(ESMF_KIND_R8), pointer         :: coordX(:,:), coordY(:,:)

    integer                             :: finalrc, result
    character(ESMF_MAXSTR)              :: testname
    character(ESMF_MAXSTR)              :: failMsg
!   !Set finalrc to success
    finalrc = ESMF_SUCCESS

    write(failMsg, *) "Example failure"
    write(testname, *) "Example ESMF_XGridEx"

    call ESMF_Initialize(defaultlogfilename="XGridEx.Log", &
                    logkindflag=ESMF_LOGKIND_MULTI, rc=localrc)
    if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT)

!>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%
!-------------------------------- Example -----------------------------
!>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%
!BOE
!\subsubsection{Create an XGrid from Grids then use it for regridding}
!\label{sec:xgrid:usage:xgrid_create}
!
! An {\tt ESMF\_XGrid} object can be created from Grids on either side
! of the exchange grid. Internally the
! weight matrices and index mapping are computed and stored in the XGrid, along
! with other necessary information for flux exchange calculation between
! any pair of model components used for the XGrid creation.
!
! In this example, we create an XGrid from overlapping Grids on
! either side of the XGrid. Then we perform a flux exchange from one side
! to the other side of the XGrid.
!
! We start by creating the Grids on both sides and associate coordinates with
! the Grids on the corner stagger. The Grids use global indexing and padding
! for coordinates on the corner stagger.
!
! For details of Grid creation and coordinate use,
! please refer to Grid class documentation: \ref{example:2DRegUniGrid}.
!EOE
!BOC
    ! First Grid on side A
    sideA(1) = ESMF_GridCreateNoPeriDim(maxIndex=(/20, 20/), &
      indexflag=ESMF_INDEX_GLOBAL, &
      gridEdgeLWidth=(/0,0/), gridEdgeUWidth=(/1,1/), &
      name='source Grid 1 on side A', rc=localrc)
!EOC
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, endflag=ESMF_END_ABORT)

!BOC
    ! Second Grid on side A
    sideA(2) = ESMF_GridCreateNoPeriDim(maxIndex=(/20, 10/), &
      indexflag=ESMF_INDEX_GLOBAL, &
      gridEdgeLWidth=(/0,0/), gridEdgeUWidth=(/1,1/), &
      name='source Grid 2 on side A', rc=localrc)
!EOC
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, endflag=ESMF_END_ABORT)

!BOC
    ! Allocate coordinates for Grid corner stagger
    do i = 1, 2
      call ESMF_GridAddCoord(sideA(i), staggerloc=ESMF_STAGGERLOC_CORNER, &
          rc=localrc)
!EOC
      if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
          endflag=ESMF_END_ABORT)
!BOC
    enddo
!EOC

!BOE
! Assign coordinate for the Grids on sideA at corner stagger.
!EOE
!BOC
    ! SideA first grid spans (0-20, 0-20) with 1.0x1.0 degree resolution
    ! X corner
    call ESMF_GridGetCoord(sideA(1), localDE=0, &
        staggerLoc=ESMF_STAGGERLOC_CORNER, coordDim=1, &
        farrayPtr=coordX, rc=localrc)
!EOC
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, endflag=ESMF_END_ABORT)
!BOC
    ! Y corner
    call ESMF_GridGetCoord(sideA(1), localDE=0, &
        staggerLoc=ESMF_STAGGERLOC_CORNER, coordDim=2, &
        farrayPtr=coordY, rc=localrc)
!EOC
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, endflag=ESMF_END_ABORT)
!BOC
    do i = lbound(coordX,1), ubound(coordX,1)
      do j = lbound(coordX, 2), ubound(coordX, 2)
        coordX(i,j) = (i-1)*1.0
        coordY(i,j) = (j-1)*1.0
      enddo
    enddo
!EOC

!BOC
    ! SideA second grid spans (14.3-24.3, 14.2-24.2) with 0.5x1.0 degree
    ! resolution X corner
    call ESMF_GridGetCoord(sideA(2), localDE=0, &
        staggerLoc=ESMF_STAGGERLOC_CORNER, coordDim=1, &
        farrayPtr=coordX, rc=localrc)
!EOC
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, endflag=ESMF_END_ABORT)
!BOC
    ! Y corner
    call ESMF_GridGetCoord(sideA(2), localDE=0, &
        staggerLoc=ESMF_STAGGERLOC_CORNER, coordDim=2, &
        farrayPtr=coordY, rc=localrc)
!EOC
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, endflag=ESMF_END_ABORT)
!BOC
    do i = lbound(coordX,1), ubound(coordX,1)
      do j = lbound(coordX, 2), ubound(coordX, 2)
        coordX(i,j) = 14.3+(i-1)*0.5
        coordY(i,j) = 14.2+(j-1)*1.0
      enddo
    enddo
!EOC

!BOE
! Create the destination grid on side B, only one Grid exists on side B. Also associate
! coordinate with the Grid:
!EOE
!BOC
    sideB(1) = ESMF_GridCreateNoPeriDim(maxIndex=(/30, 30/), &
      indexflag=ESMF_INDEX_GLOBAL, &
      gridEdgeLWidth=(/0,0/), gridEdgeUWidth=(/1,1/), &
      name='source Grid 1 on side B', rc=localrc)
!EOC
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
                endflag=ESMF_END_ABORT)
!BOC
    do i = 1, 1
      call ESMF_GridAddCoord(sideB(i), staggerloc=ESMF_STAGGERLOC_CORNER, &
          rc=localrc)
!EOC
      if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
      endflag=ESMF_END_ABORT)
!BOC
    enddo
!EOC

!BOC
    ! SideB grid spans (0-30, 0-30) with 1.0x1.0 degree resolution
    ! X corner
    call ESMF_GridGetCoord(sideB(1), localDE=0, &
        staggerLoc=ESMF_STAGGERLOC_CORNER, coordDim=1, &
        farrayPtr=coordX, rc=localrc)
!EOC
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, endflag=ESMF_END_ABORT)
!BOC
    ! Y corner
    call ESMF_GridGetCoord(sideB(1), localDE=0, &
        staggerLoc=ESMF_STAGGERLOC_CORNER, coordDim=2, &
        farrayPtr=coordY, rc=localrc)
!EOC
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, endflag=ESMF_END_ABORT)
!BOC
    do i = lbound(coordX,1), ubound(coordX,1)
      do j = lbound(coordX, 2), ubound(coordX, 2)
        coordX(i,j) = (i-1)*1.0
        coordY(i,j) = (j-1)*1.0
      enddo
    enddo
!EOC

!BOE
! Create an {\tt ESMF\_XGrid} object from the two lists of Grids on side A and B.
! In this example both Grids on side A overlaps with the Grid on side B. It's an error to have a Grid
! on either side that is spatially disjoint with the XGrid. Neither of the Grid on side A is
! identical to the Grid on side B. Calling the {\tt ESMF\_XGridCreate()} method is straightforward:
!EOE
!BOC
    xgrid = ESMF_XGridCreate(sideAGrid=sideA, sideBGrid=sideB, rc=localrc)
!EOC
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
                endflag=ESMF_END_ABORT)

!BOE
! Create an {\tt ESMF\_Field} on the XGrid:
!EOE
!BOC
    field = ESMF_FieldCreate(xgrid, typekind=ESMF_TYPEKIND_R8, &
                rc=localrc)
!EOC
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
                endflag=ESMF_END_ABORT)
!BOE
! Query the Field for its Fortran data pointer and its exclusive bounds:
!EOE
!BOC
    call ESMF_FieldGet(field, farrayPtr=xfarrayPtr, &
        exclusiveLBound=xlb, exclusiveUBound=xub, rc=localrc)
!EOC
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
                endflag=ESMF_END_ABORT)

!BOE
! Create src and dst Fields on side A and side B Grids.
!EOE
!BOC
    do i = 1, 2
        srcField(i) = ESMF_FieldCreate(sideA(i), &
                typekind=ESMF_TYPEKIND_R8, rc=localrc)
!EOC
        if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
          endflag=ESMF_END_ABORT)
!BOC
    enddo
    do i = 1, 1
        dstField(i) = ESMF_FieldCreate(sideB(i), &
                typekind=ESMF_TYPEKIND_R8, rc=localrc)
!EOC
        if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
          endflag=ESMF_END_ABORT)
!BOC
    enddo
!EOC
    dstFrac = ESMF_FieldCreate(sideB(1), &
      typekind=ESMF_TYPEKIND_R8, rc=localrc)
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
      endflag=ESMF_END_ABORT)
    dstFrac2 = ESMF_FieldCreate(sideB(1), &
      typekind=ESMF_TYPEKIND_R8, rc=localrc)
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
      endflag=ESMF_END_ABORT)


!BOE
!
! The current implementation requires that Grids used to generate the XGrid
! must not match, i.e. they are different either topologically or geometrically or both.
! In this example, the first source Grid is topologically identical to the destination
! Grid but their geometric coordinates are different.
!
! First we compute the regrid routehandles, these routehandles can be used repeatedly
! afterwards. Then we initialize the values in the Fields. Finally we execute the Regrid.
!
!EOE
!BOC
    ! Compute regrid routehandles. The routehandles can be used
    ! repeatedly afterwards.
    ! From A -> X
    do i = 1, 2
      call ESMF_FieldRegridStore(xgrid, srcField(i), field, &
        routehandle=rh_src2xgrid(i), rc = localrc)
!EOC
      if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
        endflag=ESMF_END_ABORT)
!BOC
    enddo
    ! from X -> B, retrieve the destination fraction Fields.
    do i = 1, 1
      call ESMF_FieldRegridStore(xgrid, field, dstField(i), &
        dstFracField=dstFrac, dstMergeFracField=dstFrac2, &
        routehandle=rh_xgrid2dst(i), rc = localrc)
!EOC
      if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
        endflag=ESMF_END_ABORT)
!BOC
    enddo

    ! Initialize values in the source Fields on side A
    do i = 1, 2
      call ESMF_FieldGet(srcField(i), farrayPtr=farrayPtr, rc=localrc)
!EOC
      if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
        endflag=ESMF_END_ABORT)
!BOC
      farrayPtr = i
    enddo
    ! Initialize values in the destination Field on XGrid
    xfarrayPtr = 0.0
    ! Initialize values in the destination Field on Side B
    do i = 1, 1
      call ESMF_FieldGet(dstField(i), farrayPtr=farrayPtr, rc=localrc)
!EOC
      if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
        endflag=ESMF_END_ABORT)
!BOC
      farrayPtr = 0.0
    enddo
!EOC
!BOE
! First we regrid from the Fields on side A to the Field on the XGrid:
!EOE
!BOC
    ! Execute regrid from A -> X
    do i = 1, 2
      call ESMF_FieldRegrid(srcField(i), field, &
        routehandle=rh_src2xgrid(i), &
        zeroregion=ESMF_REGION_SELECT, rc = localrc)
!EOC
      if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
        endflag=ESMF_END_ABORT)
!BOC
    enddo
!EOC

    ! xfarrayPtr should be all 1. at this point
    ! To get the surface integral of flux on XGrid, adjust by dst area

    !do i = xlb(1), xub(1)
    !    xfarrayPtr(i) = xfarrayPtr(i) * xgrid_area(i)
    !enddo

    !print *, '- after SMM from A -> X'
    !print *, xfarrayPtr ! should be xgrid_area

    !print *, '- B before SMM from X -> B'
    !print *, farrayPtr ! should be 0.
!BOE
! Next we regrid from the Field on XGrid to the destination Field on side B:
!EOE
!BOC
    ! Execute the regrid store
    do i = 1, 1
      call ESMF_FieldRegrid(field, dstField(i), &
        routehandle=rh_xgrid2dst(i), &
        rc = localrc)
!EOC
      if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
        endflag=ESMF_END_ABORT)
!BOC
    enddo
!EOC

    !print *, '- B after SMM from X -> B'
    !print *, farrayPtr ! should be 1/B_area

!BOE
! After the regridding calls, the routehandle can be released by calling the
! {\tt ESMF\_FieldRegridRelease()} method.
!EOE
!BOC
    do i = 1, 2
      call ESMF_FieldRegridRelease(routehandle=rh_src2xgrid(i), rc=localrc)
!EOC
      if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
        endflag=ESMF_END_ABORT)
!BOC
    enddo
    call ESMF_FieldRegridRelease(routehandle=rh_xgrid2dst(1), rc=localrc)
!EOC
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
      endflag=ESMF_END_ABORT)
!EOC

    call ESMF_XGridDestroy(xgrid, rc=localrc)
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
      endflag=ESMF_END_ABORT)

!BOE
! In the above example, we first set up all the required parameters to create an XGrid from user
! supplied input. Then we create Fields on the XGrid and the Grids on either side. Finally
! we use the {\tt ESMF\_FieldRegrid()} interface to perform a flux exchange from the source side
! to the destination side.
!EOE

!>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%>%
!BOE
!\subsubsection{Using XGrid in Earth System modeling}
!\label{sec:xgrid:usage:xgrid_create_masking}
!
! A typical application in Earth System Modeling is to calculate flux exchange
! through the planetary boundary layer that can be represented by {\tt ESMF\_XGrid}.
! Atmosphere is above the planetary boundary layer while land and ocean are below the boundary layer.
! To create an XGrid, the land and ocean Grids that are usually different in resolution
! need to be merged first to create a super Mesh. This merging process is enabled through the support
! of masking.
!
! The global land and ocean Grids need to be created with masking enabled.
! In practice, each Grid cell has an integer masking value attached to it. For examples using masking in
! {\tt ESMF\_Grid} please refer to section \ref{sec:usage:items}.
!
! When calling the {\tt ESMF\_XGridCreate()} method, user can supply the optional arguments
! sideAMaskValues and sideBMaskValues.
! These arguments are one dimensional Fortran integer arrays. If any of the sideAMaskValues entry
! matches the masking value used in sideA Grid, the sideA Grid cell is masked out, vice versa for sideB.
! Thus by specifying different regions of a land and ocean Grids to be masked out, the two global Grids
! can be merged into a new global Mesh covering the entire Earth.
!
! The following call shows how to use the {\tt ESMF\_XGridCreate()} method with the optional
! arguments sideAMaskValues and sideBMaskValues.
!
!EOE

!BOC
    xgrid = ESMF_XGridCreate(sideAGrid=sideA, sideBGrid=sideB, &
      sideAMaskValues=(/2/), sideBMaskValues=(/3,4/), rc=localrc)
!EOC
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
      endflag=ESMF_END_ABORT)

    call ESMF_XGridGet(xgrid, &
      sideAGridCount=ngridA, &    ! number of Grids on side A
      sideBGridCount=ngridB, &    ! number of Grids on side B
      sideAGrid=l_sideA, &    ! list of Grids on side A
      sideBGrid=l_sideB, &    ! list of Grids on side B
      distgridA=l_sideAdg, &  ! list of Distgrids on side A
      distgridM = distgrid, & ! balanced distgrid
      sparseMatA2X=l_sparseMatA2X, & !sparse matrix matmul parameters A to X
      sparseMatX2B=l_sparseMatX2B, & !sparse matrix matmul parameters X to B
      rc=localrc)
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
      endflag=ESMF_END_ABORT)

    call ESMF_XGridGet(xgrid, localDe=0, &
      elementCount=eleCount, &    ! elementCount on the localDE
      exclusiveCount=ec, &        ! exclusive count
      exclusiveLBound=elb, &      ! exclusive lower bound
      exclusiveUBound=eub, &      ! exclusive upper bound
      rc=localrc)
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
      endflag=ESMF_END_ABORT)

    call ESMF_XGridGet(xgrid, &
      xgridSide=ESMF_XGRIDSIDE_A, & ! side of the XGrid to query
      gridIndex=1, &              ! index of the distgrid
      distgrid=distgrid, &        ! the distgrid returned
      rc=localrc)
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
      endflag=ESMF_END_ABORT)
    call ESMF_XGridGet(xgrid, xgridSide=ESMF_XGRIDSIDE_A, gridIndex=2, &
      distgrid=distgrid, rc=localrc)
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
      endflag=ESMF_END_ABORT)
    call ESMF_XGridGet(xgrid, xgridSide=ESMF_XGRIDSIDE_B, gridIndex=1, &
      distgrid=distgrid, rc=localrc)
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
      endflag=ESMF_END_ABORT)

    ! After the regridding is successful.
    ! Clean up all the allocated resources:
    call ESMF_FieldDestroy(field, rc=localrc)
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
      endflag=ESMF_END_ABORT)

    call ESMF_XGridDestroy(xgrid, rc=localrc)
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
      endflag=ESMF_END_ABORT)

    do i = 1, 2
      call ESMF_FieldDestroy(srcField(i), rc = localrc)
      if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
        endflag=ESMF_END_ABORT)
      call ESMF_GridDestroy(sideA(i), rc = localrc)
      if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
        endflag=ESMF_END_ABORT)
    enddo

    do i = 1, 1
      call ESMF_FieldDestroy(dstField(i), rc = localrc)
      if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
        endflag=ESMF_END_ABORT)
      call ESMF_GridDestroy(sideB(i), rc = localrc)
      if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
        endflag=ESMF_END_ABORT)
    enddo
    call ESMF_FieldDestroy(dstFrac, rc = localrc)
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
      endflag=ESMF_END_ABORT)
    call ESMF_FieldDestroy(dstFrac2, rc = localrc)
    if(localrc /= ESMF_SUCCESS) call ESMF_Finalize(rc=localrc, &
      endflag=ESMF_END_ABORT)

    if(localrc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE
    print *, "Regridding through XGrid example returned"

    ! IMPORTANT: ESMF_STest() prints the PASS string and the # of processors in the log
    ! file that the scripts grep for.
    call ESMF_STest((finalrc.eq.ESMF_SUCCESS), testname, failMsg, result, ESMF_SRCLINE)

!-------------------------------------------------------------------------
    call ESMF_Finalize(rc=localrc)
    if (localrc.NE.ESMF_SUCCESS) finalrc = ESMF_FAILURE
!-------------------------------------------------------------------------

    if (finalrc.EQ.ESMF_SUCCESS) then
        print *, "PASS: ESMF_XGridEx.F90"
    else
        print *, "FAIL: ESMF_XGridEx.F90"
    end if
end program ESMF_XGridEx