!
!  coefh_intf1.f90
!
       subroutine coefh_intf1(cell1,cell2,ap1,ae1,aw1,an1,as1,at1,ab1,b1, &
                        h01,h1,x1,y1,z1, &
                        tc1,cp1,rho1,icvs1,jcvs1,kcvs1, &                        
                        ap2,ae2,aw2,an2,as2,at2,ab2,b2, &
                        h02,h2,x2,y2,z2, &
                        tc2,cp2,rho2,icvs2,jcvs2,kcvs2, & 
                        index1r,index2r,index1,index2,acoef1,bcoef1,anorm1,massf1, &
                        acoef2,bcoef2,anorm2,massf2,cbc,bc,ncs1,ncs2, &
                        rz,iz,dz,cz,rzp,izp,czp,ierr) 
!                       
       implicit none
       real :: rzp(*),rz(*)
       integer :: izp(*),iz(*)
       character (len=*) :: czp(*),cz(*)
       real (selected_int_kind (5)) :: dz(*)
       character (len=*) :: cell1,cell2,cbc
       real :: bc 
!
       integer :: icvs1,jcvs1,kcvs1,ierr
       real :: ap1(icvs1,jcvs1,kcvs1),ae1(icvs1,jcvs1,kcvs1),aw1(icvs1,jcvs1,kcvs1), &
               an1(icvs1,jcvs1,kcvs1),as1(icvs1,jcvs1,kcvs1), &
               at1(icvs1,jcvs1,kcvs1),ab1(icvs1,jcvs1,kcvs1),b1(icvs1,jcvs1,kcvs1), &
               h01(icvs1,jcvs1,kcvs1),h1(icvs1,jcvs1,kcvs1), &
               tc1(icvs1,jcvs1,kcvs1),cp1(icvs1,jcvs1,kcvs1),rho1(icvs1,jcvs1,kcvs1)
!
       real :: x1(icvs1,jcvs1,kcvs1),y1(icvs1,jcvs1,kcvs1),z1(icvs1,jcvs1,kcvs1)
!
       integer :: icvs2,jcvs2,kcvs2
       real :: ap2(icvs2,jcvs2,kcvs2),ae2(icvs2,jcvs2,kcvs2),aw2(icvs2,jcvs2,kcvs2), &
               an2(icvs2,jcvs2,kcvs2),as2(icvs2,jcvs2,kcvs2), &
               at2(icvs2,jcvs2,kcvs2),ab2(icvs2,jcvs2,kcvs2),b2(icvs2,jcvs2,kcvs2), &
               h02(icvs2,jcvs2,kcvs2),h2(icvs2,jcvs2,kcvs2), &
               tc2(icvs2,jcvs2,kcvs2),cp2(icvs2,jcvs2,kcvs2),rho2(icvs2,jcvs2,kcvs2)
!
       real :: x2(icvs2,jcvs2,kcvs2),y2(icvs2,jcvs2,kcvs2),z2(icvs2,jcvs2,kcvs2)
!
       integer :: ncs1,ncs2 
       real :: acoef1(ncs1),bcoef1(ncs1),anorm1(ncs1),massf1(ncs1), &
               acoef2(ncs2),bcoef2(ncs2),anorm2(ncs2),massf2(ncs2)
       integer :: index1(ncs1,4),index2(ncs2,4),index1r(6),index2r(6)
!
       real :: c1,c2,tc,bval
       real :: f,d
       real :: dtime,ap0
       integer :: i1,j1,k1,i2,j2,k2,n
       integer :: is1,js1,ks1,ie1,je1,ke1,is2,js2,ks2,ie2,je2,ke2
       logical :: least,lwest,lsouth,lnorth,lbottom,ltop,lfirst1,lfirst2                                                                                      
! ----------------------------------------------------------                    
! ---------- compute boundary coefficients for energy equation                           
       lfirst1 = .false.;lfirst2 = .false. 
       is1=index1r(1);ie1=index1r(2)
       js1=index1r(3);je1=index1r(4)
       ks1=index1r(5);ke1=index1r(6)
       is2=index2r(1);ie2=index2r(2)
       js2=index2r(3);je2=index2r(4)
       ks2=index2r(5);ke2=index2r(6)
!                                                                                         
       do n=1,ncs1                        
         i1=0; j1=0; k1=0
         if (trim(cell1).ne. 'boundary') then
           i1 = index1(n,1)
           j1 = index1(n,2)
           k1 = index1(n,3)                                                      
         endif
!
         i2=0; j2=0; k2=0           
         if (trim(cell2).ne. 'boundary') then
           i2 = index2(n,1)
           j2 = index2(n,2)
           k2 = index2(n,3)                                                                 
         endif
!
! ---------------------------------------------------------------               
! ---------- Side 1                                      
!
         if (trim(cell1).ne. 'boundary') then
!
           if ( .not.lfirst1 ) then
             lwest=.false. ;least=.false.
             if (is1==ie1) then
               if (is1==1) then
                 lwest = .true.
               else
                 least = .true.
               endif
             endif  
!
             lsouth=.false. ;lnorth=.false.
             if (js1==je1) then
               if (js1==1) then
                 lsouth = .true.
               else
                 lnorth = .true.
               endif
             endif
!
             lbottom=.false. ;ltop=.false.
             if (ks1==ke1) then
               if (ks1==1) then
                 lbottom = .true.
               else
                 ltop = .true.
               endif
             endif             
             lfirst1 = .true.
           endif
!
           f   = massf1(n)
           bval= 0.0 
!          
           if (lwest) then
             if (trim(cell2)=='boundary') then
               if (trim(cbc)=='dirichlet') then                                         
                 tc = tc1(i1,j1,k1)/cp1(i1,j1,k1)                     
                 d = tc * anorm1(n)/(x1(i1+1,j1,k1)-x1(i1,j1,k1))/2.0                          
                 bval = bc
               elseif (trim(cbc)=='zero_flux') then 
                 d = 0.0
               else
                 write(*,*) 'BC specified is unavailable'
                 stop
               endif
             else
               c1  = tc2(i2,j2,k2)/cp2(i2,j2,k2)
               c2  = tc1(i1,j1,k1)/cp1(i1,j1,k1)                        
               tc  = 2 * c1 * c2 / ( c1 + c2 )                     
               d = tc * anorm1(n)/(x1(i1,j1,k1)-x2(i2,j2,k2))
             endif
             acoef1(n) = d + f/2.0 
!
! If not a bc then set bc coef but zero it again after mg coef assembly
             if (trim(cell2)/='boundary') aw1(i1,j1,k1) = acoef1(n)
           endif
          
           if (least) then
             if (trim(cell2)=='boundary') then
               if (trim(cbc)=='dirichlet') then                                         
                 tc = tc1(i1,j1,k1)/cp1(i1,j1,k1)                     
                 d = tc * anorm1(n)/(x1(i1,j1,k1)-x1(i1-1,j1,k1))/2.0
                 bval = bc                          
               elseif (trim(cbc)=='zero_flux') then 
                 d = 0.0
               else
                 write(*,*) 'BC specified is unavailable'
                 stop
               endif           
             else
               c1 = tc2(i2,j2,k2)/cp2(i2,j2,k2)                  
               c2 = tc1(i1,j1,k1)/cp1(i1,j1,k1)                        
               tc = 2 * c1 * c2 / ( c1 + c2 )                                                                                                             
               d = tc * anorm1(n)/(x2(i2,j2,k2)-x1(i1,j1,k1))
             endif
             acoef1(n) = d + f/2.0 
!
! If not a bc then set bc coef but zero it again after mg coef assembly
             if (trim(cell2)/='boundary') ae1(i1,j1,k1) = acoef1(n)
           endif

           if (lsouth) then
             if (trim(cell2)=='boundary') then
               if (trim(cbc)=='dirichlet') then                                         
                 tc = tc1(i1,j1,k1)/cp1(i1,j1,k1)                     
                 d = tc * anorm1(n)/(y1(i1,j1+1,k1)-y1(i1,j1,k1))/2.0
                 bval = bc                          
               elseif (trim(cbc)=='zero_flux') then 
                 d = 0.0
               else
                 write(*,*) 'BC specified is unavailable'
                 stop
               endif           
             else
               c1 = tc2(i2,j2,k2)/cp2(i2,j2,k2)                 
               c2 = tc1(i1,j1,k1)/cp1(i1,j1,k1)                     
               tc = 2 * c1 * c2 / ( c1 + c2 )                     
               d = tc * anorm1(n)/(y1(i1,j1,k1)-y2(i2,j2,k2))
             endif
             acoef1(n) = d + f/2.0 
!
! If not a bc then set bc coef but zero it again after mg coef assembly
             if (trim(cell2)/='boundary') as1(i1,j1,k1) = acoef1(n)
           endif
          
           if (lnorth) then
             if (trim(cell2)=='boundary') then
               if (trim(cbc)=='dirichlet') then                                         
                 tc = tc1(i1,j1,k1)/cp1(i1,j1,k1)                     
                 d = tc * anorm1(n)/(y1(i1,j1,k1)-y1(i1,j1-1,k1))/2.0
                 bval = bc                          
               elseif (trim(cbc)=='zero_flux') then 
                 d = 0.0
               else
                 write(*,*) 'BC specified is unavailable'
                 stop
               endif           
             else
               c1 = tc2(i2,j2,k2)/cp2(i2,j2,k2)                 
               c2 = tc1(i1,j1,k1)/cp1(i1,j1,k1)                        
               tc = 2 * c1 * c2 / ( c1 + c2 )                     
               d = tc * anorm1(n)/(y2(i2,j2,k2)-y1(i1,j1,k1))
             endif
             acoef1(n) = d + f/2.0 
!
! If not a bc then set bc coef but zero it again after mg coef assembly
             if (trim(cell2)/='boundary') an1(i1,j1,k1) = acoef1(n)             
           endif
                     
           if (lbottom) then
             if (trim(cell2)=='boundary') then
               if (trim(cbc)=='dirichlet') then                                         
                 tc = tc1(i1,j1,k1)/cp1(i1,j1,k1)                     
                 d = tc * anorm1(n)/(z1(i1,j1,k1+1)-z1(i1,j1,k1))/2.0
                 bval = bc                          
               elseif (trim(cbc)=='zero_flux') then 
                 d = 0.0
               else
                 write(*,*) 'BC specified is unavailable'
                 stop
               endif            
             else
               c1 = tc2(i2,j2,k2)/cp2(i2,j2,k2)                
               c2 = tc1(i1,j1,k1)/cp1(i1,j1,k1)                        
               tc = 2 * c1 * c2 / ( c1 + c2 )           
               d = tc * anorm1(n)/(z1(i1,j1,k1)-z2(i2,j2,k2))
             endif
             acoef1(n) = d + f/2.0 
!
! If not a bc then set bc coef but zero it again after mg coef assembly
             if (trim(cell2)/='boundary') ab1(i1,j1,k1) = acoef1(n)             
           endif
          
           if (ltop) then
             if (trim(cell2)=='boundary') then
               if (trim(cbc)=='dirichlet') then                                         
                 tc = tc1(i1,j1,k1)/cp1(i1,j1,k1)                     
                 d = tc * anorm1(n)/(z1(i1,j1,k1)-z1(i1,j1,k1-1))/2.0
                 bval = bc                          
               elseif (trim(cbc)=='zero_flux') then 
                 d = 0.0
               else
                 write(*,*) 'BC specified is unavailable'
                 stop
               endif           
             else
               c1 = tc2(i2,j2,k2)/cp2(i2,j2,k2)                   
               c2 = tc1(i1,j1,k1)/cp1(i1,j1,k1)                       
               tc = 2 * c1 * c2 / ( c1 + c2 )         
               d = tc * anorm1(n)/(z2(i2,j2,k2)-z1(i1,j1,k1))
             endif
             acoef1(n) = d + f/2.0 
!
! If not a bc then set coef but zero it again after mg coef assembly
             if (trim(cell2)/='boundary') at1(i1,j1,k1) = acoef1(n)             
           endif
!
           ap1(i1,j1,k1) = ap1(i1,j1,k1) + acoef1(n)
!
           if (trim(cell2)/='boundary') then 
             b1(i1,j1,k1) = b1(i1,j1,k1) + acoef1(n)*h2(i2,j2,k2)
           else
! This assumes only heat conduction at boundary             
             b1(i1,j1,k1) = b1(i1,j1,k1) + acoef1(n)*bval*cp1(i1,j1,k1)
           endif
           bcoef1(n)= 0.0
! 
         endif
!
       enddo
!          
! ---------------------------------------------------------------               
! ---------- Side 2                                      
!
!                                                                                         
       do n=1,ncs2                        
         i1=0; j1=0; k1=0
         if (trim(cell1).ne. 'boundary') then
           i1 = index1(n,1)
           j1 = index1(n,2)
           k1 = index1(n,3)                                                      
         endif
!
         i2=0; j2=0; k2=0           
         if (trim(cell2).ne. 'boundary') then
           i2 = index2(n,1)
           j2 = index2(n,2)
           k2 = index2(n,3)                                                                 
         endif
!         
         if (trim(cell2) .ne. 'boundary') then 
!
           if ( .not.lfirst2 ) then
             lwest=.false. ;least=.false.
             if (is2==ie2) then
               if (is2==1) then
                 lwest = .true.
               else
                 least = .true.
               endif
             endif  
!
             lsouth=.false. ;lnorth=.false.
             if (js2==je2) then
               if (js2==1) then
                 lsouth = .true.
               else
                 lnorth = .true.
               endif
             endif
!
             lbottom=.false. ;ltop=.false.
             if (ks2==ke2) then
               if (ks2==1) then
                 lbottom = .true.
               else
                 ltop = .true.
               endif
             endif             
             lfirst2 = .true.
           endif
!
           f   = massf2(n)
           bval= 0.0
!          
           if (lwest) then
             if (trim(cell1)=='boundary') then
               if (trim(cbc)=='dirichlet') then                                         
                 tc = tc2(i2,j2,k2)/cp2(i2,j2,k2)                     
                 d = tc * anorm2(n)/(x2(i2+1,j2,k2)-x2(i2,j2,k2))/2.0
                 bval = bc                          
               elseif (trim(cbc)=='zero_flux') then 
                 d = 0.0
               else
                 write(*,*) 'BC specified is unavailable'
                 stop
               endif           
             else           
               c1 = tc1(i1,j1,k1)/cp1(i1,j1,k1)                 
               c2 = tc2(i2,j2,k2)/cp2(i2,j2,k2)                        
               tc = 2 * c1 * c2 / ( c1 + c2 )                     
               d = tc * anorm2(n)/(x2(i2,j2,k2)-x1(i1,j1,k1))
             endif
             acoef2(n) = d + f/2.0 
!
! If not a bc then set bc coef but zero it again after mg coef assembly
             if (trim(cell1)/='boundary') aw2(i2,j2,k2) = acoef2(n)
           endif
          
           if (least) then
             if (trim(cell1)=='boundary') then
               if (trim(cbc)=='dirichlet') then                                         
                 tc = tc2(i2,j2,k2)/cp2(i2,j2,k2)                     
                 d = tc * anorm2(n)/(x2(i2,j2,k2)-x2(i2-1,j2,k2))/2.0
                 bval = bc                          
               elseif (trim(cbc)=='zero_flux') then 
                 d = 0.0
               else
                 write(*,*) 'BC specified is unavailable'
                 stop
               endif            
             else
               c1 = tc1(i1,j1,k1)/cp1(i1,j1,k1)                  
               c2 = tc2(i2,j2,k2)/cp2(i2,j2,k2)                        
               tc = 2 * c1 * c2 / ( c1 + c2 )                                                                                                             
               d = tc * anorm2(n)/(x1(i1,j1,k1)-x2(i2,j2,k2))
             endif
             acoef2(n) = d + f/2.0 
!
! If not a bc then set bc coef but zero it again after mg coef assembly
             if (trim(cell1)/='boundary') ae2(i2,j2,k2) = acoef2(n)             
           endif

           if (lsouth) then
             if (trim(cell1)=='boundary') then
               if (trim(cbc)=='dirichlet') then                                         
                 tc = tc2(i2,j2,k2)/cp2(i2,j2,k2)                     
                 d = tc * anorm2(n)/(y2(i2,j2+1,k2)-y2(i2,j2,k2))/2.0
                 bval = bc                          
               elseif (trim(cbc)=='zero_flux') then 
                 d = 0.0
               else
                 write(*,*) 'BC specified is unavailable'
                 stop
               endif            
             else
               c1 = tc1(i1,j1,k1)/cp1(i1,j1,k1)                 
               c2 = tc2(i2,j2,k2)/cp2(i2,j2,k2)                     
               tc = 2 * c1 * c2 / ( c1 + c2 )                     
               d = tc * anorm2(n)/(y2(i2,j2,k2)-y1(i1,j1,k1))
             endif
             acoef2(n) = d + f/2.0 
!
! If not a bc then set bc coef but zero it again after mg coef assembly
             if (trim(cell1)/='boundary') as2(i2,j2,k2) = acoef2(n)             
           endif
          
           if (lnorth) then
             if (trim(cell1)=='boundary') then
               if (trim(cbc)=='dirichlet') then                                     
                 tc = tc2(i2,j2,k2)/cp2(i2,j2,k2)                     
                 d = tc * anorm2(n)/(y2(i2,j2,k2)-y2(i2,j2-1,k2))/2.0
                 bval = bc                          
               elseif (trim(cbc)=='zero_flux') then
                 d = 0.0
               else
                 write(*,*) 'BC specified is unavailable'
                 stop
               endif
             else
               c1 = tc1(i1,j1,k1)/cp1(i1,j1,k1)                 
               c2 = tc2(i2,j2,k2)/cp2(i2,j2,k2)                        
               tc = 2 * c1 * c2 / ( c1 + c2 )                     
               d = tc * anorm2(n)/(y1(i1,j1,k1)-y2(i2,j2,k2))
             endif
             acoef2(n) = d + f/2.0 
!
! If not a bc then set bc coef but zero it again after mg coef assembly
             if (trim(cell1)/='boundary') an2(i2,j2,k2) = acoef2(n)             
           endif
                     
           if (lbottom) then
             if (trim(cell1)=='boundary') then
               if (trim(cbc)=='dirichlet') then                                     
                 tc = tc2(i2,j2,k2)/cp2(i2,j2,k2)                     
                 d = tc * anorm2(n)/(z2(i2,j2,k2+1)-z2(i2,j2,k2))/2.0
                 bval = bc                          
               elseif (trim(cbc)=='zero_flux') then
                 d = 0.0
               else
                 write(*,*) 'BC specified is unavailable'
                 stop
               endif           
             else
               c1 = tc1(i1,j1,k1)/cp1(i1,j1,k1)                
               c2 = tc2(i2,j2,k2)/cp2(i2,j2,k2)                        
               tc = 2 * c1 * c2 / ( c1 + c2 )           
               d = tc * anorm2(n)/(z2(i2,j2,k2)-z1(i1,j1,k1))
             endif
             acoef2(n) = d + f/2.0 
!
! If not a bc then set bc coef but zero it again after mg coef assembly
             if (trim(cell1)/='boundary') ab2(i2,j2,k2) = acoef2(n)             
           endif
          
           if (ltop) then
             if (trim(cell1)=='boundary') then
               if (trim(cbc)=='dirichlet') then                                     
                 tc = tc2(i2,j2,k2)/cp2(i2,j2,k2)                     
                 d = tc * anorm2(n)/(z2(i2,j2,k2)-z2(i2,j2,k2-1))/2.0
                 bval = bc                          
               elseif (trim(cbc)=='zero_flux') then
                 d = 0.0
               else
                 write(*,*) 'BC specified is unavailable'
                 stop
               endif            
             else
               c1 = tc1(i1,j1,k1)/cp1(i1,j1,k1)                   
               c2 = tc2(i2,j2,k2)/cp2(i2,j2,k2)                       
               tc = 2 * c1 * c2 / ( c1 + c2 )         
               d = tc * anorm2(n)/(z1(i1,j1,k1)-z2(i2,j2,k2))
             endif
             acoef2(n) = d + f/2.0
!
! If not a bc then set coef but zero it again after mg coef assembly
             if (trim(cell1)/='boundary') at2(i2,j2,k2) = acoef2(n)
           endif 
!
           ap2(i2,j2,k2) = ap2(i2,j2,k2) + acoef2(n)
!
           if (trim(cell1)/='boundary') then 
             b2(i2,j2,k2) = b2(i2,j2,k2) + acoef2(n)*h1(i1,j1,k1)
           else  
             b2(i2,j2,k2) = b2(i2,j2,k2) + acoef2(n)*bval*cp2(i2,j2,k2)
           endif
           bcoef2(n) = 0.0
!             
         endif 
!
       enddo
                                                                                                                             
       end subroutine coefh_intf1