!
!  coefh_mg.f90
!
        subroutine coefh_mg(cell,p, &
                       ap,ae,aw,an,as,at,ab,b, & 
                       ap1,ae1,aw1,an1,as1,at1,ab1,b1, & 
                       nl,nlx,nm,nmx,nn,nnx, &                                                                                                                                                                                                                 
                       icvs,jcvs,kcvs, & 
                       icvs1,jcvs1,kcvs1, &                                              
                       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=*) :: cell 
!
        integer :: icvs,jcvs,kcvs,icvs1,jcvs1,kcvs1,p,ierr
        integer, dimension(1:10) :: nl,nlx,nm,nmx,nn,nnx
        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)
        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)
!
        integer :: l,m,n,i,j,k,is,js,ks,ie,je,ke
        real :: sumae,sumaw,suman,sumas,sumat,sumab,sumap,dp
        real :: aestar,aehat,awstar,awhat,anstar,anhat,asstar,ashat, &
                atstar,athat,abstar,abhat

!
! Initialize needed coefficients
!                     
        ae1=0.0; aw1=0.0; an1=0.0; as1=0.0               
        at1=0.0; ab1=0.0; ap1=0.0;
!
! Generate coefficients
!      
! Create and fill coarse level arrays using
!   fine grid coefficients
!
! Loop over coarse mesh
!               
        do l=1,icvs1
          do m=1,jcvs1
            do n=1,kcvs1
!
              sumae = 0.0
              sumaw = 0.0
              suman = 0.0
              sumas = 0.0
              sumat = 0.0
              sumab = 0.0
              sumap = 0.0
!
! Gather coefficients from fine mesh
!
              if (icvs == 1) then
                is = 1
                ie = icvs 
              else if (l <= nl(p)) then
                is = (l-1)*2**p+1
                ie = l*2**p
              else
                is = (l-1)*2**p+1 
                ie = is + nlx(p)-1
              endif
!
              if (jcvs == 1) then
                js = 1
                je = jcvs               
              else if (m <= nm(p)) then
                js = (m-1)*2**p+1
                je = m*2**p
              else
                js = (m-1)*2**p+1
                je = js + nmx(p)-1
              endif
!
              if (kcvs == 1) then
                ks = 1
                ke = kcvs                 
              else if (n <= nn(p)) then               
                ks = (n-1)*2**p+1
                ke = n*2**p  
              else
                ks = (n-1)*2**p+1
                ke = ks + nnx(p)-1                              
              endif
!
! Identify coefficients internal to block              
!
              do i = is,ie
              do j = js,je
              do k = ks,ke
                if (i == ie) then
                  aestar = 0.0
                  aehat  = ae(i,j,k)
                else
                  aestar = ae(i,j,k)
                  aehat  = 0.0                
                endif
                if (i == is) then
                  awstar = 0.0
                  awhat  = aw(i,j,k)
                else
                  awstar = aw(i,j,k)
                  awhat  = 0.0                
                endif
!
                if (j == je) then                 
                  anstar = 0.0
                  anhat  = an(i,j,k)                
                else
                  anstar = an(i,j,k)
                  anhat  = 0.0                 
                endif
                if (j == js) then                 
                  asstar = 0.0
                  ashat  = as(i,j,k)                
                else
                  asstar = as(i,j,k)
                  ashat  = 0.0                 
                endif
!                
                if (k == ke) then                
                  atstar = 0.0
                  athat  = at(i,j,k)
                else
                  atstar = at(i,j,k)
                  athat  = 0.0
                endif
                if (k == ks) then                
                  abstar = 0.0
                  abhat  = ab(i,j,k)
                else
                  abstar = ab(i,j,k)
                  abhat  = 0.0
                endif
!
                dp =ap(i,j,k)-aestar-awstar-anstar-asstar-atstar-abstar                   
!                                
                sumae  = sumae + aehat
                sumaw  = sumaw + awhat
                suman  = suman + anhat
                sumas  = sumas + ashat
                sumat  = sumat + athat
                sumab  = sumab + abhat
!                                                    
                sumap  = sumap + dp            
              enddo
              enddo
              enddo
!
! Apply aggregated coefficients to coarse mesh
!
              ae1(l,m,n) = sumae  
              aw1(l,m,n) = sumaw 
              an1(l,m,n) = suman  
              as1(l,m,n) = sumas               
              at1(l,m,n) = sumat  
              ab1(l,m,n) = sumab 
              ap1(l,m,n) = sumap              
!                                           
            enddo
          enddo
        enddo
!
      
        end subroutine coefh_mg