!
!  coefh.f90
!                                                                               
       subroutine coefh(cell,ap,ae,aw,an,as,at,ab, & 
                       b,h0,h,wrk,vol,anorm,x,y,z, &
                       massf,tc,cp,rho,icvs,jcvs,kcvs, &                        
                       rz,iz,dz,cz,rzp,izp,czp,ierr)                        
!                                                                                
! ------------------------------------------------------------                  
!
! subroutine to compute (structured topology):                                                        
!         - compute coefficients for energy equation                             
!          - input cv volumes
!          - face area normals
!          - volume averaged gradients (source terms)
!          - advection and diffusion schemes
! 
! ------------------------------------------------------------                  
!                                                                                
       implicit none
       real :: rzp(*),rz(*)
       integer :: izp(*),iz(*)
       character (len=*) :: czp(*),cz(*)
       real (selected_int_kind (5)) :: dz(*)
       character (len=*) :: cell 
!
       integer :: icvs,jcvs,kcvs,ierr
       real :: ap(icvs,jcvs,kcvs),ae(icvs,jcvs,kcvs),aw(icvs,jcvs,kcvs), &
               an(icvs,jcvs,kcvs),as(icvs,jcvs,kcvs),at(icvs,jcvs,kcvs), &
               ab(icvs,jcvs,kcvs),b(icvs,jcvs,kcvs),h0(icvs,jcvs,kcvs), &
               h(icvs,jcvs,kcvs),wrk(icvs,jcvs,kcvs),vol(icvs,jcvs,kcvs), &
               anorm(icvs+1,jcvs+1,kcvs+1,3),massf(icvs+1,jcvs+1,kcvs+1), &
               tc(icvs,jcvs,kcvs),cp(icvs,jcvs,kcvs),rho(icvs,jcvs,kcvs)
!
       real :: x(icvs,jcvs,kcvs),y(icvs,jcvs,kcvs),z(icvs,jcvs,kcvs)
!
       real :: tc1,tc2,tccpe,tccpw,tccpn,tccps,tccpt,tccpb
       real :: fe,fw,fn,fs,ft,fb,de,dw,dn,ds,dt,db
       real :: dtime,ap0
       integer :: i,j,k,sdiff
                                                                                        
! ----------------------------------------------------------                    
! ---------- compute coefficients for energy equation         
!
       call get_real_value('physics','/simulation/time_step', &
                           dtime,rzp,ierr) 
                                                                    
       do 100 i=1,icvs                                                      
       do 110 j=1,jcvs                                                      
       do 120 k=1,kcvs                                                      
                                                                                                                               
! ---------------------------------------------------------------               
! ---------- calculate interface flux                                     
!
          fe=0.0; fw=0.0; fn=0.0; fs=0.0; ft=0.0; fb=0.0
          de=0.0; dw=0.0; dn=0.0; ds=0.0; dt=0.0; db=0.0
!                                                                                          
          if (i < icvs) then                                                                      
            tc1   = tc(i+1,j,k)/cp(i+1,j,k)                  
            tc2   = tc(i,j,k)/cp(i,j,k)                        
            tccpe = 2 * tc1 * tc2 / ( tc1 + tc2 )
            de = tccpe * anorm(i,j,k,1)/(x(i+1,j,k) - x(i,j,k))                                 
            fe  = massf(i,j,k)
          endif
          
          if (i > 1) then                                                                      
            tc1   = tc(i-1,j,k)/cp(i-1,j,k)                 
            tc2   = tc(i,j,k)/cp(i,j,k)                        
            tccpw = 2 * tc1 * tc2 / ( tc1 + tc2 )
            dw = tccpw * anorm(i-1,j,k,1)/(x(i,j,k) - x(i-1,j,k))                                 
            fw  = massf(i-1,j,k)
          endif

          if (j < jcvs) then                                                                                    
            tc1   = tc(i,j+1,k)/cp(i,j+1,k)                 
            tc2   = tc(i,j,k)/cp(i,j,k)                        
            tccpn = 2 * tc1 * tc2 / ( tc1 + tc2 )
            dn = tccpn * anorm(i,j,k,2)/(y(i,j+1,k) - y(i,j,k))                                  
            fn  = massf(i,j,k)  
          endif
          
          if (j > 1) then                                                                      
            tc1   = tc(i,j-1,k)/cp(i,j-1,k)                 
            tc2   = tc(i,j,k)/cp(i,j,k)                     
            tccps = 2 * tc1 * tc2 / ( tc1 + tc2 )
            ds = tccps * anorm(i,j-1,k,2)/(y(i,j,k) - y(i,j-1,k))
            fs  = massf(i,j-1,k)                                 
          endif
          
          if (k < kcvs) then                                                                      
            tc1   = tc(i,j,k+1)/cp(i,j,k+1)                   
            tc2   = tc(i,j,k)/cp(i,j,k)                       
            tccpt = 2 * tc1 * tc2 / ( tc1 + tc2 )
            dt = tccpt * anorm(i,j,k,3)/(z(i,j,k+1) - z(i,j,k))
            ft  = massf(i,j,k) 
          endif                                   
          
          if (k > 1) then                                                                      
            tc1   = tc(i,j,k-1)/cp(i,j,k-1)                
            tc2   = tc(i,j,k)/cp(i,j,k)                        
            tccpb = 2 * tc1 * tc2 / ( tc1 + tc2 )
            db = tccpb * anorm(i,j,k-1,3)/(z(i,j,k) - z(i,j,k-1)) 
            fb    = massf(i,j,k-1)                                                                       
          endif  
                                                                                                                                        
! -------- energy coefficients using power law scheme                           

          ae(i,j,k) = de+fe/2.0                                                                                                        
          aw(i,j,k) = dw+fw/2.0                                                                                                                                      
          an(i,j,k) = dn+fn/2.0                                                                                                                                                                
          as(i,j,k) = ds+fs/2.0                                                                                                                    
          at(i,j,k) = dt+ft/2.0                                                                                                          
          ab(i,j,k) = db+fb/2.0 
                                                                                
          ap(i,j,k) = ae(i,j,k) + aw(i,j,k) + &                               
                      an(i,j,k) + as(i,j,k) + &                               
                      at(i,j,k) + ab(i,j,k)                                  
                                                                                                                                                                                                                                               
! ---- calculate transient term                                                 
                                      
         ap0 = vol(i,j,k) * rho(i,j,k) / dtime                                     
                                                                                                                                               
! ---- add transient term                                                       

         ap(i,j,k)  = ap(i,j,k) + ap0                                       
                                                                                
! ---- diffusion source terms                                             
                                        
         sdiff = 0.0                                                         
         b(i,j,k) = sdiff                                  
                                                                                
! ----- add transient terms                                                     
                                       
         b(i,j,k) = b(i,j,k) + ap0 * h0(i,j,k)         
                                                                                                                                                 
120    continue                                                                 
110    continue                                                                 
100    continue                                                                                                                                  
                                                                
       end subroutine coefh                                                                     
