!
!  interface_gen.f90
!
       subroutine interface_gen(is1,ie1,js1,je1,ks1,ke1, &
                                is2,ie2,js2,je2,ks2,ke2, &
				                index1,index2,index12,index21,fr12,fr21,ncs1,ncs2, &
                                rz,iz,dz,cz,rzp,izp,czp,ierr) 
!
! Testing case: interfaces are on the same plane
!               interfaces involve meshes with uniform spacing in x,y and z
!               the meshes may have different resolutions
!               interfaces share same global coordinate with i-x, j-y and k-z associations. 
!               meshes at interfaces fully overlap each other  
!
! Calculate: i) number of control surfaces (ncs)  (ncsest is a conservative estimate) 
!           ii) side 1 viewing side 2: number, node and fraction list (create cs list) 
!          iii) side 2 viewing side 1: number, node and fraction list (check cs list)
!               Check: the same number of cs should exist from either view
!
!
! Logic: i) process side 1 one face at a time
!       ii) find all faces on side 2 connected to it
!      iii) count the control surfaces created
!       iv) find the side 2 nodes and fraction of cv face   
!        v) for each cs create nodes for side1, number of nodes for side 2, nodes for side 2,
!           fraction of control-volume face for side 2 nodes
!       vi) repeat for side 2
!
       implicit none
       real :: rzp(*),rz(*)
       integer :: izp(*),iz(*),ierr
       character (len=*) :: czp(*),cz(*)
       real (selected_int_kind (5)) :: dz(*)
!
       integer :: index1(ncs1,3),index2(ncs2,3),index12(ncs1,1,3),index21(ncs2,1,3)              
       integer :: is1,ie1,js1,je1,ks1,ke1, &
                  is2,ie2,js2,je2,ks2,ke2,ncs1,ncs2
       real :: fr12(ncs1,1),fr21(ncs2,1) 
!
! Local variables
!
       integer :: i,j,k,istep1,istep2,jstep1,jstep2,kstep1,kstep2,ncschk1,ncschk2, &
                  n1,n2,n,l
!
! Tree search - future implementation
!
       istep1 = 0
       istep2 = 0 
       if (is1==ie1) then
         istep1 = is1 - is2       
         istep2 = is2 - is1
       endif	 
!
       jstep1 = 0
       jstep2 = 0 
       if (js1==je1) then
         jstep1 = js1 - js2       
         jstep2 = js2 - js1
       endif
!       	 
       kstep1 = 0
       kstep2 = 0 
       if (ks1==ke1) then
         kstep1 = ks1 - ks2       
         kstep2 = ks2 - ks1
       endif	 
!
! Start with side 1
       n = 0
       ncschk1 = 0 
       do i=is1,ie1
         do j=js1,je1
	       do k=ks1,ke1
             n = n + 1	      
!             call intersect(i,j,k,n2,list,flist)  
             n2 = 1 
             index1(n,4) = n2
	         ncschk1 = ncschk1 + n2
	         do l=1,n2
!	           index12(n,l,1) = list(l,1)  
!	           index12(n,l,2) = list(l,2)
!	           index12(n,l,3) = list(l,3)
!	           frac12(n,l) = flist(l)	       	      
               index12(n,l,1) = i + istep2
	           index12(n,l,2) = j + jstep2
	           index12(n,l,3) = k + kstep2 
	           fr12(n,l) = 1.0   
             enddo
           enddo
         enddo
       enddo 
!
! Next side 2
       n = 0
       ncschk2 = 0 
       do i=is2,ie2
         do j=js2,je2
	       do k=ks2,ke2
             n = n + 1	      
!             call intersect(i,j,k,n1,list,flist)  
             n1 = 1
             index2(n,4) = n1
	         ncschk2 = ncschk2 + n1
	         do l=1,n1
!	           index21(n,l,1) = list(l,1)  
!	           index21(n,l,2) = list(l,2)
!	           index21(n,l,3) = list(l,3)
!	           frac21(n,l) = flist(l)	       	      
	           index21(n,l,1) = i + istep1  
	           index21(n,l,2) = j + jstep1
	           index21(n,l,3) = k + kstep1
	           fr21(n,l) = 1.0
             enddo
           enddo
         enddo
       enddo
!
! Check that ncs2 = ncs1
!
       end subroutine interface_gen       
