module utilities

use mesh
use gauss

implicit none

      integer, allocatable :: ssit(:,:),its(:,:)
      double precision, allocatable :: norm(:,:)
      logical, allocatable :: bdl(:)
      double precision, allocatable :: vect(:,:),normt(:,:)
      integer, allocatable :: cip(:,:)

 contains

!
!***************************************************************
! calcul de la normal pour les faces des triangles
!***************************************************************
!
!
subroutine cal_rotnf
      implicit none
      integer i,ne1,ne2
      double precision nx,ny,d

      allocate(norm(2,nf))

      do i=1,nf
	 ne1=pseg(1,i)
	 ne2=pseg(2,i)
	 nx= (ygr(ne2)-ygr(ne1))
	 ny=-(xgr(ne2)-xgr(ne1))
	 norm(1,i)= nx
	 norm(2,i)= ny
      enddo

      return
end subroutine cal_rotnf

!
!
!******************************************************************
!
subroutine cal_its
	implicit none
	integer i,j,tk

       allocate(ssit(2,nf),its(3,ne))

        do i=1,nf
         ssit(1,i)=0
         ssit(2,i)=0
        enddo

	do i=1,nf
	   tk = tseg(1,i)
	   j = 1
	   do while (itm(j,tk).ne.i.and.j.lt.4)
	    j = j + 1
	   enddo
	   if (j.eq.4) then
	      write(*,*) 'problem in cal_its. bad connectivity table'
		write(*,*) 'face',i,' triangle',tk,' index in triangle',j
		stop
	   endif
	   its(j,tk) = 1
         ssit(1,i)=j
	   tk = tseg(2,i)
	   if (tk.gt.0) then
	   j = 1
	   do while (itm(j,tk).ne.i)
	    j = j + 1
	   enddo
	   its(j,tk) = -1
         ssit(2,i)=j
	   endif
	enddo

	return
end subroutine cal_its

!
!******************************************************************
! defines other pointers relevant to the spectral element computation
!******************************************************************
!
subroutine calpre
      implicit none

      integer nbdl
      double precision ax,bx,ay,by,deta,detai,dnx,dny,dis
      integer i,ns1,i1,i2,i3

      allocate(bdl(ne))
      allocate(vect(5,ne),normt(2,ne))

      do i=1,ne
	  bdl(i)=.false.
      enddo

      do i=1,nf
	   if (tseg(2,i).eq.0) then
	      bdl(tseg(1,i))=.true.
	   endif
      enddo

      do i=1,ne
         i1 = in(1,i)
         i2 = in(2,i)
         i3 = in(3,i)
         ax = (xgr(i2)-xgr(i1))*0.5d0
         bx = (xgr(i3)-xgr(i1))*0.5d0
         ay = (ygr(i2)-ygr(i1))*0.5d0
         by = (ygr(i3)-ygr(i1))*0.5d0
         deta = ax * by - ay * bx
         detai = 1.d0/deta
         vect(1,i)=ax
         vect(2,i)=bx
         vect(3,i)=ay
         vect(4,i)=by
         vect(5,i)=detai
         ns1 = itm(1,i)
         dnx = norm(1,ns1)
         dny = norm(2,ns1)
         dis = sqrt(dnx**2+dny**2)
         dnx = dnx / dis
         dny = dny / dis
         normt(1,i)=dnx
         normt(2,i)=dny
      enddo

      nbdl=0
      do i=1,ne
         if (bdl(i)) then
            nbdl=nbdl+1
         endif
      enddo
      write(*,*) 'nbdl:',nbdl
 
      return
end subroutine calpre

!
!
!******************************************************************
!
	subroutine coriolis(fc,f0,beta)
	implicit none
	double precision fc(nm,*),f0,beta
	integer cell,n1,n2,n3
	
	do cell=1,ne
	   n1 = in(1,cell)
	   n2 = in(2,cell)
	   n3 = in(3,cell)
	   fc(1,cell) = f0 + 0.5d0 * beta *( ygr(n3) + ygr(n2))
	   fc(2,cell) =      0.5d0 * beta * ( ygr(n3) - ygr(n1))
	   fc(nc+2,cell) =   0.5d0 * beta * ( ygr(n2) - ygr(n1))
	enddo
	return
	end

!
!*****************************************************************
! function CALIP computes the index in matrix for transformation
!*****************************************************************
!
subroutine calip

	implicit none
	integer ip,i,j,nc2

      allocate(cip(2,0:nc))

	do i=0,nc
	  ip = 0
	  j=0
	  do while(j.lt.i)
	    ip = ip + nc+1-j
	    j=j+1
	  enddo
	  cip(1,i)=ip
	enddo

	nc2=nc-1
	do i=0,nc2
	  ip = 0
	  j=0
	  do while(j.lt.i)
	    ip = ip + nc2+1-j
	    j=j+1
	  enddo
	  cip(2,i)=ip
	enddo

	return
end


end module utilities