!
!  solveint.f90
!
       subroutine solveint(p,index1r,index2r,index1,index2,index12,index21,& 
                         acoef1,bcoef1,acoef2,bcoef2,phic1,b1,phic2,b2,&
                         icvs1,jcvs1,kcvs1,icvs2,jcvs2,kcvs2,ncs1,ncs2,&
                         lbsave,ierr)
!
       implicit none
       integer :: ierr
!
       integer :: p,n,icvs1,jcvs1,kcvs1,icvs2,jcvs2,kcvs2,ncs1,ncs2
       integer :: i1,j1,k1,i2,j2,k2
       integer :: index1(ncs1,4),index2(ncs2,4),index1r(6),index2r(6),&
                  index12(ncs1,3),index21(ncs2,3)
       logical :: lbsave
!       
       real :: acoef1(ncs1),bcoef1(ncs1),acoef2(ncs2),bcoef2(ncs2)
!       
       real :: phic1(icvs1,jcvs1,kcvs1),b1(icvs1,jcvs1,kcvs1), &
               phic2(icvs1,jcvs1,kcvs1),b2(icvs1,jcvs1,kcvs1)
!
! Side 1
       if (ncs1>1 .and. ncs2 > 1) then  
         do n = 1,ncs1
!
           i1 = index1(n,1)
           j1 = index1(n,2)
           k1 = index1(n,3)
           i2 = index12(n,1)
           j2 = index12(n,2)
           k2 = index12(n,3)
!
! Adjust coefficients
           if (lbsave) bcoef1(n) = b1(i1,j1,k1)
           b1(i1,j1,k1) = bcoef1(n)+acoef1(n)*phic2(i2,j2,k2)
!           
         enddo
!
! Side 2                 
         do n = 1,ncs2
!
           i2 = index2(n,1)
           j2 = index2(n,2)
           k2 = index2(n,3)
           i1 = index21(n,1)
           j1 = index21(n,2)
           k1 = index21(n,3)
!
! Adjust coefficients
           if (lbsave) bcoef2(n) = b2(i2,j2,k2)
           b2(i2,j2,k2) = bcoef2(n)+acoef2(n)*phic1(i1,j1,k1)
!
         enddo           
       endif
!
       end subroutine solveint        