!
!  residlsc.f90
!
      subroutine residlsc(p,mi1,mj1,mk1, &                           
                         ap1,ae1,aw1,an1,as1,at1,ab1,&
                         b1,phi1,mi2,mj2,mk2,r2o,res)                                        
!                                                                                                 
      implicit none
!
      integer :: p,mi1,mj1,mk1
      real :: ap1(mi1,mj1,mk1),ae1(mi1,mj1,mk1),aw1(mi1,mj1,mk1), &
              an1(mi1,mj1,mk1),as1(mi1,mj1,mk1),                  &
              at1(mi1,mj1,mk1),ab1(mi1,mj1,mk1)
      real :: b1(mi1,mj1,mk1),phi1(mi1,mj1,mk1)
      integer :: mi2,mj2,mk2 
      real :: r2o(mi2,mj2,mk2)
!
      integer :: isum,lend1,mend1,nend1 
      integer :: i,j,k,l,m,n,is,ie,js,je,ks,ke
      real :: res,sumr               
!
      real :: sumae,sumaw,suman,sumas,sumat,sumab,sumap
      real :: aestar,awstar,anstar,asstar,atstar,abstar
      real :: aehat,awhat,anhat,ashat,athat,abhat,dp 
!                                                                                      
! ------------------------------------------------------------                  
!
! Smooth residuals on a coarser mesh
!
      res = 0.0   
      isum= 0 
!
      if (p>=1) then
!
! Loop over coarse mesh (p points to level of coarse mesh)
!
        lend1 = mod(mi1,2)
        mend1 = mod(mj1,2) 
        nend1 = mod(mk1,2)       
!              
        do l=1,mi2
          do m=1,mj2
            do n=1,mk2
!                 
              sumr  = 0.0
!
! Gather range for the fine mesh block
!
              if (mi1 == 1) then
                is = 1
                ie = mi1 
              else if (l < mi2) then
                is = (l-1)*2+1
                ie = is+1
              else
                is = (l-1)*2+1
                ie = is+1 - lend1
              endif
!
              if (mj1 == 1) then
                js = 1
                je = mj1               
              else if (m < mj2) then
                js = (m-1)*2+1
                je = js + 1
              else
                js = (m-1)*2+1
                je = js+1 - mend1
              endif
!
              if (mk1 == 1) then
                ks = 1
                ke = mk1                 
              else if (n < mk2) then               
                ks = (n-1)*2+1
                ke = ks + 1  
              else
                ks = (n-1)*2+1
                ke = ks+1 - nend1                                 
              endif
!
! Sum residuals from fine mesh block              
!
              do i = is,ie
              do j = js,je
              do k = ks,ke
!                                    
                if (i == ie) then
                  aestar = 0.0
                  aehat  = ae1(i,j,k)
                else
                  aestar = ae1(i,j,k)
                  aehat  = 0.0                
                endif
                if (i == is) then
                  awstar = 0.0
                  awhat  = aw1(i,j,k)
                else
                  awstar = aw1(i,j,k)
                  awhat  = 0.0                
                endif
!
                if (j == je) then                 
                  anstar = 0.0
                  anhat  = an1(i,j,k)                
                else
                  anstar = an1(i,j,k)
                  anhat  = 0.0                 
                endif
                if (j == js) then                 
                  asstar = 0.0
                  ashat  = as1(i,j,k)                
                else
                  asstar = as1(i,j,k)
                  ashat  = 0.0                 
                endif
!                
                if (k == ke) then                
                  atstar = 0.0
                  athat  = at1(i,j,k)
                else
                  atstar = at1(i,j,k)
                  athat  = 0.0
                endif
                if (k == ks) then                
                  abstar = 0.0
                  abhat  = ab1(i,j,k)
                else
                  abstar = ab1(i,j,k)
                  abhat  = 0.0
                endif
!
                dp=ap1(i,j,k)-aestar-awstar-anstar-asstar-atstar-abstar                   
!                                
                if (i==mi1) then
                  sumae = 0.0
                else 
                  sumae  = aehat*phi1(i+1,j,k)
                endif
                if (i==1) then  
                  sumaw  = 0.0
                else
                  sumaw  = awhat*phi1(i-1,j,k)
                endif
!                  
                if (j==mj1) then
                  suman  = 0.0
                else
                  suman  = anhat*phi1(i,j+1,k)
                endif
                if (j==1) then  
                  sumas = 0.0
                else
                  sumas  = ashat*phi1(i,j-1,k)
                endif
!
                if (k==mk1) then
                  sumat = 0.0 
                else                   
                  sumat = athat*phi1(i,j,k+1)
                endif
                if (k==1) then  
                  sumab  = 0.0
                else
                  sumab  = abhat*phi1(i,j,k-1) 
                endif 
!                                                    
                sumap  = dp*phi1(i,j,k) 
!
                sumr=sumr+b1(i,j,k)+ &
                     sumae+sumaw+suman+sumas+sumat+sumab-sumap 
!
              enddo
              enddo
              enddo
!
! Apply aggregated residuals to coarse mesh
!
              r2o(l,m,n) = sumr               
              res   = res  + abs(sumr)
              isum  = isum + 1
!                                           
            enddo
          enddo
        enddo      
!
      else
!
        write(*,*) 'Incorrect level in residlsc'
        stop
!          
      endif 
!
      res = res/isum
!                                                                                                                                                                 
      end subroutine residlsc