!     Compute geometric arrays used in SELFE (side info etc)
!     In the driver routine, the following 2 routines should be called
!     in the following fashion (all float arguments are real*4):
!
!     First get # of sides 'ns' with inputs: np (# of nodes),ne (# of elements), and nm (connectivity table)
!     call compute_nside(np,ne,nm,ns) 
!     Allocate side-related arrays
!     allocate(ic3(ne,3),js(ne,3),is(ns,2),isidenode(ns,2),xcj(ns,2),ycj(ns,2)
!     Then compute the rest of side related arrays with additional inputs (xnd,ynd) (x,y coordinates of each node)
!     call selfe_geometry(np,ne,ns,xnd,ynd,nm,ic3,js,is,isidenode,xcj,ycj)
!     The outputs:
!                ns: # of sides
!                ic3(ne,3): 3 neighboring elements of an element
!                js(ne,3): 3 sides of an element
!                is(ns,2): 2 adjacent elements of a side
!                isidenode(ns,2): 2 end nodes of a side
!                xcj(ns),ycj(ns): x,y of each side center

!     Common data used in 2 routines
      module selfe_geometry_mod
      implicit none
      public
      integer,save :: nx(4,4,3)
      integer,save,allocatable :: nne(:),ine(:,:),ic3(:,:)
      end module selfe_geometry_mod

      subroutine compute_nside(np,ne,nm,ns)
      use selfe_geometry_mod
      implicit real(4)(a-h,o-z),integer(i-n)
      integer, intent(in) :: np,ne,nm(ne,3)
      integer, intent(out) :: ns !,ic3(ne,3)
!      integer, allocatable :: ine(:,:)

      allocate(nne(np),ic3(ne,3),stat=istat)
      if(istat/=0) stop 'Failed to alloc. nne'

!     Compute geometry
      do k=3,4
        do i=1,k
          do j=1,k-1
            nx(k,i,j)=i+j
            if(nx(k,i,j)>k) nx(k,i,j)=nx(k,i,j)-k
            if(nx(k,i,j)<1.or.nx(k,i,j)>k) then
              write(*,*)'nx wrong',i,j,k,nx(k,i,j)
              stop
            endif
          enddo !j
        enddo !i
      enddo !k

      nne=0
      do i=1,ne
        do j=1,3
          nd=nm(i,j)
          nne(nd)=nne(nd)+1
!          ine(nd,nne(nd))=i
        enddo
      enddo
      mnei=maxval(nne)

      allocate(ine(np,mnei),stat=istat)
      if(istat/=0) stop 'Failed to alloc. ine'
      nne=0
      do i=1,ne
        do j=1,3
          nd=nm(i,j)
          nne(nd)=nne(nd)+1
          if(nne(nd)>mnei) then
            write(*,*)'Too many neighbors',nd
            stop
          endif
          ine(nd,nne(nd))=i
        enddo
      enddo !i

!     Compute ball info; this won't be affected by re-arrangement below
      do i=1,ne
        do j=1,3 
          ic3(i,j)=0 !index for bnd sides
          nd1=nm(i,nx(3,j,1))
          nd2=nm(i,nx(3,j,2))
          do k=1,nne(nd1)
            ie=ine(nd1,k)
            if(ie/=i.and.(nm(ie,1)==nd2.or.nm(ie,2)==nd2.or.nm(ie,3)==nd2)) ic3(i,j)=ie
          enddo !k
        enddo !j
      enddo !i

      ns=0
      do ie=1,ne
        do j=1,3 !visit each side associated with element ie
          if(ic3(ie,j)==0.or.ie<ic3(ie,j)) then !new side
            ns=ns+1
          endif
        enddo !j
      enddo !ie

      end subroutine compute_nside

      subroutine selfe_geometry(np,ne,ns0,xnd,ynd,nm,ic3_out,js,is,isidenode,xcj,ycj)
      use selfe_geometry_mod
      implicit real(4)(a-h,o-z),integer(i-n)
      integer, intent(in) :: np,ne,ns0,nm(ne,3)
      real, intent(in) :: xnd(np),ynd(np)
      integer, intent(out) :: ic3_out(ne,3),js(ne,3),is(ns0,2),isidenode(ns0,2)
      real, intent(out) :: xcj(ns0),ycj(ns0)
      
      ic3_out=ic3

      ns=0 !# of sides
      do i=1,ne
        do j=1,3
          nd1=nm(i,nx(3,j,1))
          nd2=nm(i,nx(3,j,2))
          if(ic3(i,j)==0.or.i<ic3(i,j)) then !new sides
            ns=ns+1
            if(ns>ns0) then
              write(*,*)'Too many sides'
              stop
            endif
            js(i,j)=ns
            is(ns,1)=i
            isidenode(ns,1)=nd1
            isidenode(ns,2)=nd2
            xcj(ns)=(xnd(nd1)+xnd(nd2))/2
            ycj(ns)=(ynd(nd1)+ynd(nd2))/2

            is(ns,2)=ic3(i,j) !bnd element => bnd side
!           Corresponding side in element ic3(i,j)
            if(ic3(i,j)/=0) then !old internal side
              iel=ic3(i,j)
              index=0
              do k=1,3
                if(ic3(iel,k)==i) then
                  index=k
                  exit
                endif
              enddo !k
              if(index==0) then
                write(*,*)'Wrong ball info',i,j
                stop
              endif
              js(iel,index)=ns
            endif !ic3(i,j).ne.0
          endif !ic3(i,j)==0.or.i<ic3(i,j)
        enddo !j
      enddo !i=1,ne

      if(ns/=ns0) stop 'Side count mismatch'

      end subroutine selfe_geometry