!
!  solvegs_seq.f90
!
      subroutine solvegs_seq(ap,ae,aw,an,as,at,ab,b,phi,mi,mj,mk)                    
!                                                                                
      implicit none                                                                                                                       
!                                                                                
      real :: phi(mi,mj,mk), ap(mi,mj,mk), ae(mi,mj,mk), &                     
              aw(mi,mj,mk), an(mi,mj,mk), as(mi,mj,mk),  &                    
              at(mi,mj,mk), ab(mi,mj,mk), b(mi,mj,mk)
!                      
      integer :: mi,mj,mk
!
      real :: aephi,awphi,anphi,asphi,atphi,abphi
      integer :: i,j,k
!                                                                              
      do i=1, mi
        do j=1, mj  
          do k=1, mk 
!
            if (i==mi) then
              aephi = 0.0
            else                                                        
              aephi = ae(i,j,k)*phi(i+1,j,k)
            endif
            if (i==1) then
              awphi = 0.0
            else                                                        
              awphi = aw(i,j,k)*phi(i-1,j,k)
            endif          
!
            if (j==mj) then
              anphi = 0.0
            else                                                        
              anphi = an(i,j,k)*phi(i,j+1,k)
            endif
            if (j==1) then
              asphi = 0.0
            else                                                        
              asphi = as(i,j,k)*phi(i,j-1,k)
            endif 
!                                                                                                 
            if (k==mk) then
              atphi = 0.0
            else                                                        
              atphi = at(i,j,k)*phi(i,j,k+1)
            endif
            if (k==1) then
              abphi = 0.0
            else                                                        
              abphi = ab(i,j,k)*phi(i,j,k-1) 
            endif 
!                                                                                                              
            phi(i,j,k)=                &
              (aephi + awphi + anphi + asphi + atphi + abphi +  &
               b(i,j,k))/ap(i,j,k)
          enddo                                                             
        enddo                                                                
      enddo                                                                  
!                                                                                                                                                    
      end subroutine solvegs_seq 
