 !
!  interface_input.f90
!
       subroutine interface_input(rz,iz,dz,cz,rzp,izp,czp,ierr)      
!
       implicit none
       real :: rzp(*),rz(*)
       integer :: izp(*),iz(*),ierr
       character (len=*) :: czp(*),cz(*)
       real (selected_int_kind (5)) :: dz(*)
       character (len=80) :: cval
!
       integer :: n,m,nint,ntmp,ncs1,ncs2,np
       integer :: icvs,jcvs,kcvs,icvp,jcvp,kcvp,length
       integer :: i,j,k,p,is1,ie1,js1,je1,ks1,ke1,is2,ie2,js2,je2,ks2,ke2
       character (len=80) :: intf,cell1,cell2,intp,cell1p,cell2p 
       character (len=4) :: ctype
       character (len=2), dimension(1:10) :: cnum
       data cnum/'1 ','2 ','3 ','4 ','5 ','6 ','7 ','8 ','9 ','10'/
       logical :: least1,lwest1,lnorth1,lsouth1,ltop1,lbottom1, &
                  least2,lwest2,lnorth2,lsouth2,ltop2,lbottom2 
!
! Pointers
       integer :: pint,pcind,pindr1,pindr2,ptrr,ptrc,pcbc,pbc
       integer :: pan1,pan2,pmf1,pmf2,ptr1,ptr2,ptr12,ptr21,pfr12,pfr21
       integer :: px,py,pz
!
       call get_mem('physics','/interface_list',ctype,pint,nint,ierr)
!
       do n=1,nint
         intf = '/'//czp(pint + n -1) 
         call get_mem('physics',trim(intf)//'/cell_index',ctype,pcind,ntmp,ierr)
         cell1=czp(pcind)
         cell2=czp(pcind+1)
!
! Determine cell ranges: at present only structured cells
         if (trim(cell1) .ne. 'boundary') then
           call get_mem('physics',trim(intf)//'/index1_range',ctype,pindr1,ntmp,ierr)                 
           is1 = izp(pindr1)
           ie1 = izp(pindr1+1)
           js1 = izp(pindr1+2)
           je1 = izp(pindr1+3)         
           ks1 = izp(pindr1+4)
           ke1 = izp(pindr1+5)
           ncs1= (ie1-is1+1)*(je1-js1+1)*(ke1-ks1+1) 
         else
           is1 = 1; ie1 = 1; js1 = 1; je1 = 1; ks1 = 1; ke1 = 1
           ncs1= 1         
         endif
!
         if (trim(cell2) .ne. 'boundary') then
           call get_mem('physics',trim(intf)//'/index2_range',ctype,pindr2,ntmp,ierr)         
           is2 = izp(pindr2)
           ie2 = izp(pindr2+1)
           js2 = izp(pindr2+2)
           je2 = izp(pindr2+3)
           ks2 = izp(pindr2+4)
           ke2 = izp(pindr2+5)
           ncs2= (ie2-is2+1)*(je2-js2+1)*(ke2-ks2+1) 
         else
           is2 = 1; ie2 = 1; js2 = 1; je2 = 1; ks2 = 1; ke2 = 1
           ncs2= 1
         endif             
!
         call get_next_ptr('interface','int',ptr1,ierr)
         call put_mem('interface',trim(intf)//'/index1','int',ptr1,ncs1*4,ierr)
         call mkintf_list(iz(ptr1),ncs1,is1,ie1,js1,je1,ks1,ke1)        
!
         call get_next_ptr('interface','int',ptr2,ierr)
         call put_mem('interface',trim(intf)//'/index2','int',ptr2,ncs2*4,ierr)
         call mkintf_list(iz(ptr2),ncs2,is2,ie2,js2,je2,ks2,ke2)
!
         if (trim(cell1)=='boundary' .or. trim(cell2)=='boundary') then       
!         
           call get_next_ptr('interface','real',ptrr,ierr)         
           call put_mem('interface',trim(intf)//'/acoef1','real',ptrr,ncs1,ierr)
           call set_a_0(rz(ptrr),ncs1)
           ptrr = ptrr + ncs1 
           call put_mem('interface',trim(intf)//'/bcoef1','real',ptrr,ncs1,ierr)
           call set_a_0(rz(ptrr),ncs1)
           pan1 = ptrr + ncs1 
           call put_mem('interface',trim(intf)//'/anorm1','real',pan1,ncs1,ierr)
           call set_a_0(rz(pan1),ncs1)
           pmf1 = pan1 + ncs1 
           call put_mem('interface',trim(intf)//'/massfl1','real',pmf1,ncs1,ierr)
           call set_a_0(rz(pmf1),ncs1)
!         
           call get_next_ptr('interface','real',ptrr,ierr)         
           call put_mem('interface',trim(intf)//'/acoef2','real',ptrr,ncs2,ierr) 
           call set_a_0(rz(ptrr),ncs2)
           ptrr = ptrr + ncs2 
           call put_mem('interface',trim(intf)//'/bcoef2','real',ptrr,ncs2,ierr)
           call set_a_0(rz(ptrr),ncs2)
           pan2 = ptrr + ncs2 
           call put_mem('interface',trim(intf)//'/anorm2','real',pan2,ncs2,ierr)
           call set_a_0(rz(pan2),ncs2)
           pmf2 = pan2 + ncs2 
           call put_mem('interface',trim(intf)//'/massfl2','real',pmf2,ncs2,ierr)
           call set_a_0(rz(pmf2),ncs2)
!
! Fill in control surface areas and mass flows
           if (trim(cell1) .ne. 'boundary') then                   
             call get_int_value('physics',trim(cell1)//'/icvs',icvs,izp,ierr)  
             call get_int_value('physics',trim(cell1)//'/jcvs',jcvs,izp,ierr)      
             call get_int_value('physics',trim(cell1)//'/kcvs',kcvs,izp,ierr)
!
             call get_mem('cell',trim(cell1)//'/mesh/x',ctype,px,length,ierr)
             call get_mem('cell',trim(cell1)//'/mesh/y',ctype,py,length,ierr)
             call get_mem('cell',trim(cell1)//'/mesh/z',ctype,pz,length,ierr)
!                                    
             call cal_anorm_intf(is1,ie1,js1,je1,ks1,ke1,iz(ptr1), &
                                 rz(px),rz(py),rz(pz),rz(pan1),&
                                 icvs+1,jcvs+1,kcvs+1,ncs1,ierr)
! Walls with zero mass flow for now
           else              
             call get_int_value('physics',trim(cell2)//'/icvs',icvs,izp,ierr)  
             call get_int_value('physics',trim(cell2)//'/jcvs',jcvs,izp,ierr)      
             call get_int_value('physics',trim(cell2)//'/kcvs',kcvs,izp,ierr)
!
             call get_mem('cell',trim(cell2)//'/mesh/x',ctype,px,length,ierr)
             call get_mem('cell',trim(cell2)//'/mesh/y',ctype,py,length,ierr)
             call get_mem('cell',trim(cell2)//'/mesh/z',ctype,pz,length,ierr)                                    
!
             call cal_anorm_intf(is2,ie2,js2,je2,ks2,ke2,iz(ptr2), &
                                 rz(px),rz(py),rz(pz),rz(pan2),&
                                 icvs+1,jcvs+1,kcvs+1,ncs2,ierr)               
! Walls with zero mass flow for now
           endif           
         else  
!
! If there is no boundary then find intersection between interfaces
!
           call get_next_ptr('interface','int',ptr12,ierr)
           call put_mem('interface',trim(intf)//'/index12','int',ptr12,ncs1*3,ierr)       
!
           call get_next_ptr('interface','int',ptr21,ierr)
           call put_mem('interface',trim(intf)//'/index21','int',ptr21,ncs2*3,ierr)
!
           call get_next_ptr('interface','real',pfr12,ierr)
           call put_mem('interface',trim(intf)//'/frac12','int',pfr12,ncs1*3,ierr)       
!
           call get_next_ptr('interface','real',pfr21,ierr)
           call put_mem('interface',trim(intf)//'/frac21','int',pfr21,ncs2*3,ierr)
! 
           call interface_gen(is1,ie1,js1,je1,ks1,ke1, &
                              is2,ie2,js2,je2,ks2,ke2, &
		                      iz(ptr1),iz(ptr2),iz(ptr12),iz(ptr21),rz(pfr12),rz(pfr21), &
                              ncs1,ncs2,rz,iz,dz,cz,rzp,izp,czp,ierr) 
!
           call get_next_ptr('interface','real',ptrr,ierr)         
           call put_mem('interface',trim(intf)//'/acoef1','real',ptrr,ncs1,ierr) 
           call set_a_0(rz(ptrr),ncs1)
           ptrr = ptrr + ncs1 
           call put_mem('interface',trim(intf)//'/bcoef1','real',ptrr,ncs1,ierr)
           call set_a_0(rz(ptrr),ncs1)
           pan1 = ptrr + ncs1 
           call put_mem('interface',trim(intf)//'/anorm1','real',pan1,ncs1,ierr)
           call set_a_0(rz(pan1),ncs1)
           pmf1 = pan1 + ncs1 
           call put_mem('interface',trim(intf)//'/massfl1','real',pmf1,ncs1,ierr)
           call set_a_0(rz(pmf1),ncs1)
!           
           call get_int_value('physics',trim(cell1)//'/icvs',icvs,izp,ierr)  
           call get_int_value('physics',trim(cell1)//'/jcvs',jcvs,izp,ierr)      
           call get_int_value('physics',trim(cell1)//'/kcvs',kcvs,izp,ierr)
!
           call get_mem('cell',trim(cell1)//'/mesh/x',ctype,px,length,ierr)
           call get_mem('cell',trim(cell1)//'/mesh/y',ctype,py,length,ierr)
           call get_mem('cell',trim(cell1)//'/mesh/z',ctype,pz,length,ierr)
!                                    
           call cal_anorm_intf(is1,ie1,js1,je1,ks1,ke1,iz(ptr1), &
                               rz(px),rz(py),rz(pz),rz(pan1),&
                               icvs+1,jcvs+1,kcvs+1,ncs1,ierr)
!
! Define space for ip mass flows
!                                
           call get_next_ptr('interface','real',ptrr,ierr)         
           call put_mem('interface',trim(intf)//'/acoef2','real',ptrr,ncs2,ierr) 
           call set_a_0(rz(ptrr),ncs2)
           ptrr = ptrr + ncs2 
           call put_mem('interface',trim(intf)//'/bcoef2','real',ptrr,ncs2,ierr)
           call set_a_0(rz(ptrr),ncs2)
           pan2 = ptrr + ncs2 
           call put_mem('interface',trim(intf)//'/anorm2','real',pan2,ncs2,ierr)
           call set_a_0(rz(pan2),ncs2)
           pmf2 = pan2 + ncs2 
           call put_mem('interface',trim(intf)//'/massfl2','real',pmf2,ncs2,ierr)
           call set_a_0(rz(pmf2),ncs2)
!           
           call get_int_value('physics',trim(cell2)//'/icvs',icvs,izp,ierr)  
           call get_int_value('physics',trim(cell2)//'/jcvs',jcvs,izp,ierr)      
           call get_int_value('physics',trim(cell2)//'/kcvs',kcvs,izp,ierr)
!
           call get_mem('cell',trim(cell2)//'/mesh/x',ctype,px,length,ierr)
           call get_mem('cell',trim(cell2)//'/mesh/y',ctype,py,length,ierr)
           call get_mem('cell',trim(cell2)//'/mesh/z',ctype,pz,length,ierr)                                    
!
           call cal_anorm_intf(is2,ie2,js2,je2,ks2,ke2,iz(ptr2), &
                               rz(px),rz(py),rz(pz),rz(pan2),&
                               icvs+1,jcvs+1,kcvs+1,ncs2,ierr) 
!
! Define space for ip mass flows
!
         endif          
       enddo  
!
! Create multigrid level interfaces
!
! --Loop over interfaces
!   - Grab the cell names on each side of the interface
!
! --Loop over the mg level interfaces
!   - Check if one side is a boundary condition
!   - For each cell get the icvs, jcvs and kcvs values 
!   - Create the int_#/level_#/index1_range etc. (as done in setup_physics)
!   - If no bc involved then calculate the number of control surfaces, else
!     number of control surfaces is set by the one adjacent mesh.
!   - Create data areas 
!   - Repeat
!
       do n=1,nint
         intf = '/'//czp(pint + n -1) 
         call get_mem('physics',trim(intf)//'/cell_index',ctype,pcind,ntmp,ierr)
         cell1=czp(pcind)
         cell2=czp(pcind+1)
!
! Determine cell ranges: at present only structured cells
         least1=.false.;lwest1=.false.;lnorth1=.false.;lsouth1=.false.;ltop1=.false.;lbottom1=.false. 
         if (trim(cell1) .ne. 'boundary') then
           call get_mem('physics',trim(intf)//'/index1_range',ctype,pindr1,ntmp,ierr)                 
           is1 = izp(pindr1)
           ie1 = izp(pindr1+1)
           if (is1==ie1) then
             if (is1==1) then 
               lwest1 = .true.
             else
               least1 = .true.
             endif  
           endif 
!
           js1 = izp(pindr1+2)
           je1 = izp(pindr1+3)         
           if (js1==je1) then
             if (js1==1) then 
               lsouth1 = .true.
             else
               lnorth1 = .true.
             endif  
           endif 
!
           ks1 = izp(pindr1+4)
           ke1 = izp(pindr1+5)
           if (ks1==ke1) then
             if (ks1==1) then 
               lbottom1 = .true.
             else
               ltop1 = .true.
             endif  
           endif                     
         endif
!
         least2=.false.;lwest2=.false.;lnorth2=.false.;lsouth2=.false.;ltop2=.false.;lbottom2=.false. 
         if (trim(cell2) .ne. 'boundary') then
           call get_mem('physics',trim(intf)//'/index2_range',ctype,pindr2,ntmp,ierr)
           is2 = izp(pindr2)
           ie2 = izp(pindr2+1)
           if (is2==ie2) then
             if (is2==1) then 
               lwest2 = .true.
             else
               least2 = .true.
             endif  
           endif 
!
           js2 = izp(pindr2+2)
           je2 = izp(pindr2+3)
           if (js2==je2) then
             if (js2==1) then 
               lsouth2 = .true.
             else
               lnorth2 = .true.
             endif  
           endif 
!
           ks2 = izp(pindr2+4)
           ke2 = izp(pindr2+5)
           if (ks2==ke2) then
             if (ks2==1) then 
               lbottom2 = .true.
             else
               ltop2 = .true.
             endif  
           endif            
         endif
!
         call get_int_value('physics','/simulation/npcells',np,izp,ierr)
         do p=1,np
           intp  = trim(intf)//'/level_'//cnum(p) 
!
           if (trim(cell1) .ne. 'boundary') then 
             cell1p= trim(cell1)//'/level_'//cnum(p)
           else 
             cell1p= cell1
           endif  
           if (trim(cell2) .ne. 'boundary') then
             cell2p= trim(cell2)//'/level_'//cnum(p)
           else
             cell2p= cell2
           endif
!
! Determine cell ranges: at present only structured cells
           if (trim(cell1) .ne. 'boundary') then
             call get_int_value('coef',trim(cell1p)//'/icvs',icvp,iz,ierr)
             call get_int_value('coef',trim(cell1p)//'/jcvs',jcvp,iz,ierr)
             call get_int_value('coef',trim(cell1p)//'/kcvs',kcvp,iz,ierr) 	 
             is1 = 1; ie1 = icvp
             if (lwest1) ie1 = 1; if (least1) is1 = ie1              
             js1 = 1; je1 = jcvp
             if (lsouth1) je1 = 1; if (lnorth1) js1 = je1                      
             ks1 = 1; ke1 = kcvp
             if (lbottom1) ke1 = 1; if (ltop1) ks1 = ke1 
!
             call get_next_ptr('interface','int',pindr1,ierr)	 
             call put_mem('interface',trim(intp)//'/index1_range','int',pindr1,6,ierr)                 
             iz(pindr1)  = is1
             iz(pindr1+1)= ie1
             iz(pindr1+2)= js1
             iz(pindr1+3)= je1         
             iz(pindr1+4)= ks1
             iz(pindr1+5)= ke1
             ncs1= (ie1-is1+1)*(je1-js1+1)*(ke1-ks1+1) 
           else
             is1 = 1; ie1 = 1; js1 = 1; je1 = 1; ks1 = 1; ke1 = 1
             ncs1= 1         
           endif
!
           if (trim(cell2) .ne. 'boundary') then
             call get_int_value('coef',trim(cell2p)//'/icvs',icvp,iz,ierr)
             call get_int_value('coef',trim(cell2p)//'/jcvs',jcvp,iz,ierr)
             call get_int_value('coef',trim(cell2p)//'/kcvs',kcvp,iz,ierr) 
             is2 = 1; ie2 = icvp
             if (lwest2) ie2 = 1; if (least2) is2 = ie2              
             js2 = 1; je2 = jcvp
             if (lsouth2) je1 = 1; if (lnorth2) js2 = je2                      
             ks2 = 1; ke2 = kcvp
             if (lbottom2) ke2 = 1; if (ltop2) ks2 = ke2 
!
             call get_next_ptr('interface','int',pindr2,ierr)	   
             call put_mem('interface',trim(intp)//'/index2_range','int',pindr2,6,ierr)         
             iz(pindr2)  = is2
             iz(pindr2+1)= ie2
             iz(pindr2+2)= js2
             iz(pindr2+3)= je2
             iz(pindr2+4)= ks2
             iz(pindr2+5)= ke2
             ncs2= (ie2-is2+1)*(je2-js2+1)*(ke2-ks2+1) 
           else
             is2 = 1; ie2 = 1; js2 = 1; je2 = 1; ks2 = 1; ke2 = 1
             ncs2= 1
           endif             
!
           call get_next_ptr('interface','int',ptr1,ierr)
           call put_mem('interface',trim(intp)//'/index1','int',ptr1,ncs1*4,ierr)
           call mkintf_list(iz(ptr1),ncs1,is1,ie1,js1,je1,ks1,ke1)        
!
           call get_next_ptr('interface','int',ptr2,ierr)
           call put_mem('interface',trim(intp)//'/index2','int',ptr2,ncs2*4,ierr)
           call mkintf_list(iz(ptr2),ncs2,is2,ie2,js2,je2,ks2,ke2)
!         
           call get_next_ptr('interface','real',ptrr,ierr)         
           call put_mem('interface',trim(intp)//'/acoef1','real',ptrr,ncs1,ierr) 
           call set_a_0(rz(ptrr),ncs1)
           ptrr = ptrr + ncs1 
           call put_mem('interface',trim(intp)//'/bcoef1','real',ptrr,ncs1,ierr)
           call set_a_0(rz(ptrr),ncs1)
!         
           call get_next_ptr('interface','real',ptrr,ierr)         
           call put_mem('interface',trim(intp)//'/acoef2','real',ptrr,ncs2,ierr) 
           call set_a_0(rz(ptrr),ncs2)
           ptrr = ptrr + ncs2 
           call put_mem('interface',trim(intp)//'/bcoef2','real',ptrr,ncs2,ierr)
           call set_a_0(rz(ptrr),ncs2)
!
           if (trim(cell1)/='boundary' .and. trim(cell2)/='boundary') then 
!
! If there is no boundary then find intersection between interfaces
!
             call get_next_ptr('interface','int',ptr12,ierr)
             call put_mem('interface',trim(intp)//'/index12','int',ptr12,ncs1*3,ierr)       
!
             call get_next_ptr('interface','int',ptr21,ierr)
             call put_mem('interface',trim(intp)//'/index21','int',ptr21,ncs2*3,ierr)
!
             call get_next_ptr('interface','real',pfr12,ierr)
             call put_mem('interface',trim(intp)//'/frac12','int',pfr12,ncs1*3,ierr)       
!
             call get_next_ptr('interface','real',pfr21,ierr)
             call put_mem('interface',trim(intp)//'/frac21','int',pfr21,ncs2*3,ierr)
! 
             call interface_gen(is1,ie1,js1,je1,ks1,ke1, &
                                is2,ie2,js2,je2,ks2,ke2, &
				                iz(ptr1),iz(ptr2),iz(ptr12),iz(ptr21),rz(pfr12),rz(pfr21), &
                                ncs1,ncs2,rz,iz,dz,cz,rzp,izp,czp,ierr)             
!
           endif          
         enddo 
       enddo
!         
       end subroutine interface_input                   
