calculateConnect Subroutine

subroutine calculateConnect(minIndexPTile, maxIndexPTile, tilepair, &
           contactTuple, orientationVector, positionVector)

integer, pointer :: minIndexPTile(:,:)
integer, pointer :: maxIndexPTile(:,:)
integer, intent(IN) :: tilepair(2)
integer, intent(IN) :: contactTuple(2,4)
integer, intent(OUT) :: orientationVector(2)
integer, intent(OUT) :: positionVector(2)


  integer     :: minIndexPair(2,2)
  integer     :: maxIndexPair(2,2)
  integer     :: pA(2)   ! tile A start point
  integer     :: pB(2)   ! tile A end point
  integer     :: pC(2)   ! tile B start point (connect with point A)
  integer     :: pD(2)   ! tile B end point (connect with point B)
  ! -- helper variables for conversion
  integer     :: i, dimCount=2
  integer     :: a(2), b(2)
  integer     :: pAplus(2), pAplusrot(2)
  integer     :: a_zero, a_notzero
  integer     :: b_zero, b_notzero

  minIndexPair(:,1)=minIndexPTile(:,tilepair(1))
  minIndexPair(:,2)=minIndexPTile(:,tilepair(2))
  maxIndexPair(:,1)=maxIndexPTile(:,tilepair(1))
  maxIndexPair(:,2)=maxIndexPTile(:,tilepair(2))

  ! -- resulting DistGrid orientation vector and position vector aruments
  ! --- convesion from pA, pB, pC, pD  --> orientationVector, positionVector
  pA = contactTuple(:,1)
  pB = contactTuple(:,2)
  pC = contactTuple(:,3)
  pD = contactTuple(:,4)

  a = pB - pA ! contact vector on tile A
  b = pD - pC ! contact vector on tile B
  !print *, "a=", a
  !print *, "b=", b

  a_zero    = -1 ! invalidate
  a_notzero = -1 ! invalidate
  do i=1, dimCount
    if (a(i)==0) then
      a_zero=i
    else
      a_notzero = i
    endif
  enddo
  if (a_zero==-1 .or. a_notzero==-1) then
    call ESMF_LogSetError(ESMF_RC_ARG_BAD, &
      msg="Contact for tile A must line up with index space.")
    return
  endif

  b_zero    = -1 ! invalidate
  b_notzero = -1 ! invalidate
  do i=1, dimCount
    if (b(i)==0) then
      b_zero=i
    else
      b_notzero = i
    endif
  enddo
  if (b_zero==-1 .or. b_notzero==-1) then
    call ESMF_LogSetError(ESMF_RC_ARG_BAD, &
      msg="Contact for tile B must line up with index space.")
    return
  endif

  !print *, "a_zero=", a_zero, "   a_notzero=", a_notzero
  !print *, "b_zero=", b_zero, "   b_notzero=", b_notzero

  if (abs(a(a_notzero)) /= abs(b(b_notzero))) then
    call ESMF_LogSetError(ESMF_RC_ARG_BAD, &
      msg="Contact must have same number of elements on both sides.")
    return
  endif

  ! construct preliminary orientationVector
  orientationVector(b_zero) = a_zero
  orientationVector(b_notzero) = a_notzero

  ! consider sign flip along contact
  if (a(a_notzero) == -b(b_notzero)) then
    orientationVector(b_notzero) = - orientationVector(b_notzero)
  endif

  ! consider sign flip perpendicular to contact
  if ((pA(a_zero)==minindexPair(a_zero,1) .and. &
       pC(b_zero)==minindexPair(b_zero,2)) .or. &
      (pA(a_zero)==maxIndexPair(a_zero,1) .and. &
       pC(b_zero)==maxIndexPair(b_zero,2))) then
    orientationVector(b_zero) = - orientationVector(b_zero)
  endif

  ! done constructing orientationVector
  ! print *, ">>>>> orientationVector=", orientationVector, " <<<<<"

  ! now construct positionVector....
  ! The contact connects pA "+1" (on tile A) to pC (on tile B). The "+1" here
  ! means that one step beyond tile A is taken. This step is perpendicular to
  ! the contact vector, and outward from the tile.
  pAplus = pA
  if (pA(a_zero)==minindexPair(a_zero,1)) pAplus(a_zero) = pAplus(a_zero) - 1
  if (pA(a_zero)==maxIndexPair(a_zero,1)) pAplus(a_zero) = pAplus(a_zero) + 1
  !print *, "pAplus=", pAplus
  ! now apply the previously determined orientationVector on pAplus
  do i=1, dimCount
    if (orientationVector(i) < 0) then
      pAplusrot(i) = -pAplus(-orientationVector(i))
    else
      pAplusrot(i) = pAplus(orientationVector(i))
    endif
  enddo
  !print *, "pAplusrot=", pAplusrot
  ! find positionVector per definition as the difference
  positionVector = pC - pAplusrot
  !print *, ">>>>> positionVector=   ", positionVector, " <<<<<"

end subroutine calculateConnect