ESMF_XGridCreateFromSparseMat Function

function ESMF_XGridCreateFromSparseMat(keywordEnforcer, &
    sideAGrid,              sideAMesh, &
    sideBGrid,              sideBMesh, &
    sideAGridPriority,      sideAMeshPriority, &
    sideBGridPriority,      sideBMeshPriority, &
    sparseMatA2X, sparseMatX2A, sparseMatB2X, sparseMatX2B, &
    area, centroid, &
    name, &
    rc) 

!
! !RETURN VALUE:
    type(ESMF_XGrid) :: ESMF_XGridCreateFromSparseMat
!
! !ARGUMENTS:
type(ESMF_KeywordEnforcer), optional:: keywordEnforcer ! must use keywords below
type(ESMF_Grid),      intent(in), optional :: sideAGrid(:)
type(ESMF_Mesh),      intent(in), optional :: sideAMesh(:)
type(ESMF_Grid),      intent(in), optional :: sideBGrid(:)
type(ESMF_Mesh),      intent(in), optional :: sideBMesh(:)
integer,              intent(in), optional :: sideAGridPriority(:)
integer,              intent(in), optional :: sideAMeshPriority(:)
integer,              intent(in), optional :: sideBGridPriority(:)
integer,              intent(in), optional :: sideBMeshPriority(:)
type(ESMF_XGridSpec), intent(in), optional :: sparseMatA2X(:)
type(ESMF_XGridSpec), intent(in), optional :: sparseMatX2A(:)
type(ESMF_XGridSpec), intent(in), optional :: sparseMatB2X(:)
type(ESMF_XGridSpec), intent(in), optional :: sparseMatX2B(:)
real(ESMF_KIND_R8),   intent(in), optional :: area(:)
real(ESMF_KIND_R8),   intent(in), optional :: centroid(:,:)
character (len=*),    intent(in), optional :: name
integer,              intent(out),optional :: rc 

!
! !DESCRIPTION:
!      Create an XGrid directly from user supplied sparse matrix parameters. User
!      is responsible to supply all information necessary for communication calculation. 
!      For an example and associated documentation using this method see section 
!      \ref{sec:xgrid:usage:xgrid_createfromsparsemat}
!
!     The arguments are:
!     \begin{description}
!     \item [{[sideAGrid]}]
!           Parametric 2D Grids on side A, for example, 
!           these Grids can be either Cartesian 2D or Spherical.
!     \item [{[sideAMesh]}]
!           Parametric 2D Meshes on side A, for example, 
!           these Meshes can be either Cartesian 2D or Spherical.
!     \item [{[sideBGrid]}]
!           Parametric 2D Grids on side B, for example, 
!           these Grids can be either Cartesian 2D or Spherical.
!     \item [{[sideBMesh]}]
!           Parametric 2D Meshes on side B, for example, 
!           these Meshes can be either Cartesian 2D or Spherical.
!     \item [{[sideAGridPriority]}]
!           Priority array of Grids on sideA during overlay generation.
!           The {\tt sideAGridPriority} array should be the same size as the {\tt sideAGrid} array. The values
!           in the array should range from 1 to size(sideAGrid)+size(sideAMesh). A Grid whose corresponding
!           value in this array is lower than another side A Grid or Mesh, will take precedence over that Grid or Mesh
!           during side A merging. In other words, if both have parts in the same region, then the object with the lower value will win, and
!           the other Grid or Mesh part will be clipped away.
!     \item [{[sideAMeshPriority]}]
!           Priority array of Meshes on sideA during overlay generation.
!           The {\tt sideAMeshPriority} array should be the same size as the {\tt sideAMesh} array. The values
!           in the array should range from 1 to size(sideAGrid)+size(sideAMesh). A Mesh whose corresponding
!           value in this array is lower than another side A Grid or Mesh, will take precedence over that Grid or Mesh
!           during side A merging. In other words, if both have parts in the same region, then the object with the lower value will win, and
!           the other Grid or Mesh part will be clipped away.
!     \item [{[sideBGridPriority]}]
!           Priority array of Grids on sideB during overlay generation.
!           The {\tt sideBGridPriority} array should be the same size as the {\tt sideBGrid} array. The values
!           in the array should range from 1 to size(sideBGrid)+size(sideBMesh). A Grid whose corresponding
!           value in this array is lower than another side B Grid or Mesh, will take precedence over that Grid or Mesh
!           during side B merging. In other words, if both have parts in the same region, then the object with the lower value will win, and
!           the other Grid or Mesh part will be clipped away.
!     \item [{[sideBMeshPriority]}]
!           Priority array of Meshes on sideB during overlay generation.
!           The {\tt sideBMeshPriority} array should be the same size as the {\tt sideBMesh} array. The values
!           in the array should range from 1 to size(sideBGrid)+size(sideBMesh). A Mesh whose corresponding
!           value in this array is lower than another side B Grid or Mesh, will take precedence over that Grid or Mesh
!           during side B merging. In other words, if both have parts in the same region, then the object with the lower value will win, and
!           the other Grid or Mesh part will be clipped away.
!     \item [{[sparseMatA2X]}]
!           indexlist from a Grid index space on side A to xgrid index space;
!           indexFactorlist from a Grid index space on side A to xgrid index space.
!     \item [{[sparseMatX2A]}]
!           indexlist from xgrid index space to a Grid index space on side A;
!           indexFactorlist from xgrid index space to a Grid index space on side A.
!     \item [{[sparseMatB2X]}]
!           indexlist from a Grid index space on side B to xgrid index space;
!           indexFactorlist from a Grid index space on side B to xgrid index space.
!     \item [{[sparseMatX2B]}]
!           indexlist from xgrid index space to a Grid index space on side B;
!           indexFactorlist from xgrid index space to a Grid index space on side B.
!     \item [{[area]}]
!           area of the xgrid cells.
!     \item [{[centroid]}]
!           coordinates at the area weighted center of the xgrid cells.
!     \item [{[name]}]
!           name of the xgrid object.
!     \item [{[rc]}]
!           Return code; equals {\tt ESMF\_SUCCESS} only if the {\tt ESMF\_XGrid} 
!           is created.
!     \end{description}
!
!EOP

    integer :: localrc, ngrid_a, ngrid_b
    integer :: i, j, ncells, ndim
    type(ESMF_XGridType), pointer :: xgtype

    ! Initialize
    localrc = ESMF_RC_NOT_IMPL

    ! Initialize return code   
    if(present(rc)) rc = ESMF_RC_NOT_IMPL

    ! check there are enough input to create the XGrid
    if(.not. present(sideAGrid) .and. .not. present(sideAMesh)) then
      call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, & 
         msg="- Either Grid or Mesh must be provided on sideA", &
         ESMF_CONTEXT, rcToReturn=rc) 
      return
    endif
    if(.not. present(sideBGrid) .and. .not. present(sideBMesh)) then
      call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, & 
         msg="- Either Grid or Mesh must be provided on sideB", &
         ESMF_CONTEXT, rcToReturn=rc) 
      return
    endif

    if(present(sideAGrid) .and. present(sideAMesh)) then
      if(size(sideAGrid, 1)+size(sideAMesh, 1) .le. 0) then
        call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, & 
           msg="- Either Grid or Mesh must be provided on sideA", &
           ESMF_CONTEXT, rcToReturn=rc) 
        return
      endif
    endif

    if(present(sideBGrid) .and. present(sideBMesh)) then
      if(size(sideBGrid, 1)+size(sideBMesh, 1) .le. 0) then
        call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, & 
           msg="- Either Grid or Mesh must be provided on sideB", &
           ESMF_CONTEXT, rcToReturn=rc) 
        return
      endif
    endif

    ! check init status of input Grids
    if(present(sideAGrid)) then
      ngrid_a = size(sideAGrid, 1)
      do i = 1, ngrid_a
          ESMF_INIT_CHECK_DEEP(ESMF_GridGetInit,sideAGrid(i),rc)
      enddo

      ! No need to check grid for offline xgrid creation, assume user know what they are doing!
      !do i = 1, ngrid_a
      !  call checkGrid(sideAGrid(i), ESMF_STAGGERLOC_CORNER, rc=localrc)
      !  if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
      !      ESMF_CONTEXT, rcToReturn=rc)) return
      !enddo

      if(present(sideAGridPriority)) then
        if(size(sideAGridPriority, 1) /= ngrid_a) then
          call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, & 
             msg="- Size of sideAGridPriority does not match size of sideAGrid", &
             ESMF_CONTEXT, rcToReturn=rc) 
          return
        endif
      endif
    endif

    if(present(sideBGrid)) then
      ngrid_b = size(sideBGrid, 1)
      do i = 1, ngrid_b
          ESMF_INIT_CHECK_DEEP(ESMF_GridGetInit,sideBGrid(i),rc)
      enddo
      ! No need to check grid for offline xgrid creation, assume user know what they are doing!
      !do i = 1, ngrid_b
      !  call checkGrid(sideBGrid(i), ESMF_STAGGERLOC_CORNER, rc=localrc)
      !  if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
      !      ESMF_CONTEXT, rcToReturn=rc)) return
      !enddo
      if(present(sideBGridPriority)) then
        if(size(sideBGridPriority, 1) /= ngrid_a) then
          call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, & 
             msg="- Size of sideBGridPriority does not match size of sideBGrid", &
             ESMF_CONTEXT, rcToReturn=rc) 
          return
        endif
      endif
    endif

    !TODO: need to check Meshes too.
    if(present(sideAMesh)) then
      ngrid_a = size(sideAMesh, 1)
      if(present(sideAMeshPriority)) then
        if(size(sideAMeshPriority, 1) /= ngrid_a) then
          call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, & 
             msg="- Size of sideAMeshPriority does not match size of sideAMesh", &
             ESMF_CONTEXT, rcToReturn=rc) 
          return
        endif
      endif
    endif
    if(present(sideBMesh)) then
      ngrid_b = size(sideBMesh, 1)
      if(present(sideBMeshPriority)) then
        if(size(sideBMeshPriority, 1) /= ngrid_b) then
          call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, & 
             msg="- Size of sideBMeshPriority does not match size of sideBMesh", &
             ESMF_CONTEXT, rcToReturn=rc) 
          return
        endif
      endif
    endif

    ! Priority for Grid and Mesh must all be present simultaneously
    if( (present(sideAGridPriority) .and. .not. present(sideAGrid)) .or. &
         present(sideAMeshPriority) .and. .not. present(sideAMesh)  .or. &
         present(sideBGridPriority) .and. .not. present(sideBGrid)  .or. &
         present(sideBMeshPriority) .and. .not. present(sideBMesh)) then

      call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, & 
         msg="- Cannot specify Grid or Mesh Priority without actual list of Grids or Meshes", &
         ESMF_CONTEXT, rcToReturn=rc) 
      return
    endif

    ! Make sure that the input objects have a consisent coordSys
    call error_check_coordsys(sideAGrid, sideAMesh, sideBGrid, sideBMesh, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
         ESMF_ERR_PASSTHRU, &
         ESMF_CONTEXT, rcToReturn=rc)) return


    ! At this point the inputs are consistently sized.
    ! Take care of ordering
    ngrid_a = 0
    if(present(sideAGrid)) ngrid_a = size(sideAGrid, 1)
    if(present(sideAMesh)) ngrid_a = ngrid_a + size(sideAMesh, 1)
    ngrid_b = 0
    if(present(sideBGrid)) ngrid_b = size(sideBGrid, 1)
    if(present(sideBMesh)) ngrid_b = ngrid_b + size(sideBMesh, 1)

    ! do some range checking on priority lists
    if(present(sideAGridPriority) .and. present(sideAMeshPriority)) then
      do i = 1, size(sideAGridPriority, 1)
        if(sideAGridPriority(i) .le. 0 .or. sideAGridPriority(i) .gt. ngrid_a) then
          call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, & 
             msg="- sideAGridPriority value out of range", &
             ESMF_CONTEXT, rcToReturn=rc) 
          return
        endif
      enddo
      do i = 1, size(sideAMeshPriority, 1)
        if(sideAMeshPriority(i) .le. 0 .or. sideAMeshPriority(i) .gt. ngrid_a) then
          call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, & 
             msg="- sideAMeshPriority value out of range", &
             ESMF_CONTEXT, rcToReturn=rc) 
          return
        endif
      enddo

      do i = 1, size(sideAGridPriority, 1)
        do j = 1, size(sideAMeshPriority, 1)
          if(sideAGridPriority(i) == sideAMeshPriority(j)) then
            call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, & 
               msg="- sideAGridPriority and sideAMeshPriority cannot have duplicate entry", &
               ESMF_CONTEXT, rcToReturn=rc) 
            return
          endif
        enddo
      enddo
    endif

    if(present(sideBGridPriority) .and. present(sideBMeshPriority)) then
      do i = 1, size(sideBGridPriority, 1)
        if(sideBGridPriority(i) .le. 0 .or. sideBGridPriority(i) .gt. ngrid_b) then
          call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, & 
             msg="- sideBGridPriority value out of range", &
             ESMF_CONTEXT, rcToReturn=rc) 
          return
        endif
      enddo
      do i = 1, size(sideBMeshPriority, 1)
        if(sideBMeshPriority(i) .le. 0 .or. sideBMeshPriority(i) .gt. ngrid_b) then
          call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, & 
             msg="- sideBMeshPriority value out of range", &
             ESMF_CONTEXT, rcToReturn=rc) 
          return
        endif
      enddo

      do i = 1, size(sideBGridPriority, 1)
        do j = 1, size(sideBMeshPriority, 1)
          if(sideBGridPriority(i) == sideBMeshPriority(j)) then
            call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, & 
               msg="- sideBGridPriority and sideBMeshPriority cannot have duplicate entry", &
               ESMF_CONTEXT, rcToReturn=rc) 
            return
          endif
        enddo
      enddo
    endif

    if(.not. present(sparseMatA2X) .and. &
      (.not. present(sparseMatX2A)).and. &
      (.not. present(sparseMatB2X)).and. &
      (.not. present(sparseMatX2B))) then
        call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, & 
           msg="- One of the sparseMat must be set to generate an XGrid", &
           ESMF_CONTEXT, rcToReturn=rc) 
        return
    endif

    ! initialize XGridType object and its base object
    nullify(xgtype)
    nullify(ESMF_XGridCreateFromSparseMat%xgtypep)
    call ESMF_XGridConstructBaseObj(xgtype, name, localrc)
    if (ESMF_LogFoundAllocError(localrc, &
      msg="Constructing xgtype base object ", &
      ESMF_CONTEXT, rcToReturn=rc)) return

    ! Choose a default coordSys for XGrid
    ! TODO: eventually also allow a user to pass one in
    call choose_default_coordsys(sideAGrid, sideAMesh, sideBGrid, sideBMesh, xgtype%coordSys, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
         ESMF_ERR_PASSTHRU, &
         ESMF_CONTEXT, rcToReturn=rc)) return

    allocate(xgtype%sideA(ngrid_a), xgtype%sideB(ngrid_b), stat=localrc)
    if (ESMF_LogFoundAllocError(localrc, &
        msg="- Allocating xgtype%sideA or xgtype%sideB ", &
        ESMF_CONTEXT, rcToReturn=rc)) return

    ! Need to initialize xgtype%sideA based on sideAGrid and/or sideAMesh
    if(present(sideAGrid)) then
      if(present(sideAGridPriority)) then
        do i = 1, size(sideAGrid, 1)
          xgtype%sideA(sideAGridPriority(i)) = ESMF_XGridGeomBaseCreate(sideAGrid(i), &
            ESMF_STAGGERLOC_CENTER, rc=localrc)
          if (ESMF_LogFoundError(localrc, &
              ESMF_ERR_PASSTHRU, &
              ESMF_CONTEXT, rcToReturn=rc)) return
        enddo
      else
        do i = 1, size(sideAGrid, 1)
          xgtype%sideA(i) = ESMF_XGridGeomBaseCreate(sideAGrid(i), &
            ESMF_STAGGERLOC_CENTER, rc=localrc)
          if (ESMF_LogFoundError(localrc, &
              ESMF_ERR_PASSTHRU, &
              ESMF_CONTEXT, rcToReturn=rc)) return
        enddo
      endif

      if(present(sideAMesh)) then
        if(present(sideAMeshPriority)) then
          do i = 1, size(sideAMesh, 1)
            xgtype%sideA(sideAMeshPriority(i)) = ESMF_XGridGeomBaseCreate(sideAMesh(i), &
              ESMF_MESHLOC_ELEMENT, rc=localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
          enddo
        else
          do i = 1, size(sideAMesh, 1)
            xgtype%sideA(i+size(sideAGrid, 1)) = ESMF_XGridGeomBaseCreate(sideAMesh(i), &
              ESMF_MESHLOC_ELEMENT, rc=localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
          enddo
        endif
      endif
    else ! .not. present(sideAGrid)
      if(present(sideAMesh)) then
        if(present(sideAMeshPriority)) then
          do i = 1, size(sideAMesh, 1)
            xgtype%sideA(sideAMeshPriority(i)) = ESMF_XGridGeomBaseCreate(sideAMesh(i), &
              ESMF_MESHLOC_ELEMENT, rc=localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
          enddo
        else
          do i = 1, size(sideAMesh, 1)
            xgtype%sideA(i) = ESMF_XGridGeomBaseCreate(sideAMesh(i), &
              ESMF_MESHLOC_ELEMENT, rc=localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
          enddo
        endif
      endif
    endif

    ! Need to initialize xgtype%sideB based on sideBGrid and/or sideBMesh
    if(present(sideBGrid)) then
      if(present(sideBGridPriority)) then
        do i = 1, size(sideBGrid, 1)
          xgtype%sideB(sideBGridPriority(i)) = ESMF_XGridGeomBaseCreate(sideBGrid(i), &
            ESMF_STAGGERLOC_CENTER, rc=localrc)
          if (ESMF_LogFoundError(localrc, &
              ESMF_ERR_PASSTHRU, &
              ESMF_CONTEXT, rcToReturn=rc)) return
        enddo
      else
        do i = 1, size(sideBGrid, 1)
          xgtype%sideB(i) = ESMF_XGridGeomBaseCreate(sideBGrid(i), &
            ESMF_STAGGERLOC_CENTER, rc=localrc)
          if (ESMF_LogFoundError(localrc, &
              ESMF_ERR_PASSTHRU, &
              ESMF_CONTEXT, rcToReturn=rc)) return
        enddo
      endif

      if(present(sideBMesh)) then
        if(present(sideBMeshPriority)) then
          do i = 1, size(sideBMesh, 1)
            xgtype%sideB(sideBMeshPriority(i)) = ESMF_XGridGeomBaseCreate(sideBMesh(i), &
              ESMF_MESHLOC_ELEMENT, rc=localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
          enddo
        else
          do i = 1, size(sideBMesh, 1)
            xgtype%sideB(i+size(sideBGrid, 1)) = ESMF_XGridGeomBaseCreate(sideBMesh(i), &
              ESMF_MESHLOC_ELEMENT, rc=localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
          enddo
        endif
      endif
    else ! .not. present(sideBGrid)
      if(present(sideBMesh)) then
        if(present(sideBMeshPriority)) then
          do i = 1, size(sideBMesh, 1)
            xgtype%sideB(sideBMeshPriority(i)) = ESMF_XGridGeomBaseCreate(sideBMesh(i), &
              ESMF_MESHLOC_ELEMENT, rc=localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
          enddo
        else
          do i = 1, size(sideBMesh, 1)
            xgtype%sideB(i) = ESMF_XGridGeomBaseCreate(sideBMesh(i), &
              ESMF_MESHLOC_ELEMENT, rc=localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
          enddo
        endif
      endif
    endif

    if(present(area)) then
      ncells = size(area, 1)
      allocate(xgtype%area(ncells), stat=localrc)
      if (ESMF_LogFoundAllocError(localrc, &
          msg="- Allocating xgtype%area ", &
          ESMF_CONTEXT, rcToReturn=rc)) return
      xgtype%area = area
    endif

    if(present(centroid)) then
        ncells = size(centroid, 1)
        ndim = size(centroid, 2)
        allocate(xgtype%centroid(ncells, ndim), stat=localrc)
        if (ESMF_LogFoundAllocError(localrc, &
            msg="- Allocating xgtype%centroid ", &
            ESMF_CONTEXT, rcToReturn=rc)) return
        xgtype%centroid = centroid
    endif

    call ESMF_XGridConstruct(xgtype, xgtype%sideA, xgtype%sideB, &
      sparseMatA2X=sparseMatA2X, sparseMatX2A=sparseMatX2A, &
      sparseMatB2X=sparseMatB2X, sparseMatX2B=sparseMatX2B, offline=.true., rc=localrc)
    if (ESMF_LogFoundAllocError(localrc, &
      msg="Constructing xgtype object ", &
      ESMF_CONTEXT, rcToReturn=rc)) return

    ! Finalize XGrid Creation
    xgtype%online = 0
    xgtype%status = ESMF_STATUS_READY
    ESMF_XGridCreateFromSparseMat%xgtypep => xgtype 
    ESMF_INIT_SET_CREATED(ESMF_XGridCreateFromSparseMat)

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

    if(present(rc)) rc = ESMF_SUCCESS

end function ESMF_XGridCreateFromSparseMat