c***********************************************************************
c     This subroutine invert the pressure equation:
c
c     div (1/rho grad p) = rhs
c***
c     Li's ALGEBRAIC VERSION
c***  
c     ON INPUT:
c     
c     3D arrays (nx,ny,nz):
c
c       rhs: initial right hand side (divergence field)
c       u: pressure
c     res: memory reserved for residual on the top grid
c       a: coefficients for equation at each level (computed in acfromac)
c       b: coefficients "   "
c       c: coefficients "   "
c      ei: inverse coefficients  "  "
c   
c     scalar parameters:
c 
c     nx,ny,nz: dimensions of arrays
c       cycles: number of cycles at each V level
c
c     ON OUTPUT:
c     
c     everything the same,
c***********************************************************************
      SUBROUTINE mglin()
c***
      INCLUDE 'undefined.h'
      INCLUDE 'dimension.h'
      INCLUDE 'pointers.h'
      INCLUDE 'physics.h'
      INCLUDE 'algorithms_controllers.h'
      INTEGER j,jcycle,jj,n,jpost,jpre,nf,ngrid,nn,nxmin,nymin,nzmin
      INTEGER nnx,nny,nnz,nfx,nfy,nfz
      DOUBLE PRECISION norm, conv_rate
      DOUBLE PRECISION FINDN2S
      EXTERNAL FINDN2S
      INTRINSIC MIN0,MAX0
      SAVE conv_rate
      DATA conv_rate /0.5d0/

c     SUBROUTINES: copy, rstrct, abcei, fill0, slvsml, interp,
c     relax, resid, addint
c***
c     get minimum dimension
c***
      nn  = MIN0(nx,ny,nz)
      nnx = nx
      nny = ny
      nnz = nz
      nfx = nx
      nfy = ny
      nfz = nz
      ngrid = NG

c***
c     calculate matrix coefficients on every grid
c***
      do while (nn .gt. 4)
         nn  =  nn/2 + 1
         nnx = nnx/2 + 1
         nny = nny/2 + 1
         nnz = nnz/2 + 1
         ngrid = ngrid - 1

         call acfromac(za(ia(ngrid)),zb(ia(ngrid)),zc(ia(ngrid)),
     %        zei(ia(ngrid)),za(ia(ngrid+1)),zb(ia(ngrid+1)),
     %        zc(ia(ngrid+1)),nnx,nny,nnz)
      enddo
c***
c     define minimum values, just for clarity
c***
      nn = 4
      nxmin = nnx
      nymin = nny
      nzmin = nnz

c      WRITE(*,*)cycles,NPRE,NPOST,NTOP,OMEGA

c Test MG : start from zero pressure
c
c      call fill0(zu(ia(NG)),nfx,nfy,nfz)

c***
c     Copy previous solution at time t^n in a work field d(NG)
c***
      call copy(zd(ia(NG)),zu(ia(NG)),nfx,nfy,nfz)

c***  
c     Calculate the residual on the finest grid and put it in rhs(NG) so that 
c     we will actually relax on the error eq.  Ae=res.
c***
      
      call resid(zres(ia(NG)),zd(ia(NG)),zrhs(ia(NG)),
     %         za(ia(NG)),zb(ia(NG)),zc(ia(NG)),nfx,nfy,nfz)

      call copy(zrhs(ia(NG)),zres(ia(NG)),nfx,nfy,nfz)

c***
c     if bad convergence rate at previous time step, then increase the number of 
c       relaxations for the new cycle
c***
      RELAX_TIMES_OLD=RELAX_TIMES
      if (conv_rate .gt. 0.6) then 
         RELAX_TIMES = RELAX_TIMES + 1
         write(*,*) "convergence rate deteriorating, nbr iterations = ",relax_times
       else if (conv_rate .lt. 0.1) then
         RELAX_TIMES = MAX0(RELAX_TIMES - 1,1)
         write(*,*) "convergence rate improving,     nbr iterations = ",relax_times
      endif

c Need norm to compute norm old

      norm=h**2*FINDN2S(zrhs(ia(NG)),nfx,nfy,nfz)
      if(printres.eq.1) then
         write(*,*) "mglin: begin norm =",norm
      endif

      if(NCYCLE.lt.2) STOP 'NCYCLE too small'
      do n=1,NCYCLE

c Test if we can exit loop
c
         if(norm.lt.mineps.and.n.gt.1) then
            call copy(zu(ia(NG)),zd(ia(NG)),nfx,nfy,nfz)

c  commented out by Denis ????

            call bc_scalar(zu(ia(NG)),nfx,nfy,nfz)
            conv_rate=norm/NORM_OLD
            NORM_OLD=norm
            if(printres.eq.1) then
               write(*,*) 'mglin: return cycles=',n-1,' norm=',norm
            endif  
            cycles = n-1
            return
         endif


c     on finest grid relax on the error e(NG) with initial guess e(NG)=0 

      call fill0(zu(ia(NG)),nfx,nfy,nfz)

c     Now u contains the error (or correction) and not the solution anymore.

      do jj=NG,2,-1
c***
c     compute   relax_times relaxations on the finest grid, 
c             2*relax_times on the NG-1th  grid..
c      . (NG-1)*relax_times on the coarsest grid
c***
            NPRE=RELAX_TIMES*(NG-jj+1)

            do jpre=1,NPRE
               call relax(zu(ia(jj)),zrhs(ia(jj)),za(ia(jj)),
     %              zb(ia(jj)),zc(ia(jj)),zei(ia(jj)),
     %              nfx,nfy,nfz,OMEGA)
            enddo

            call resid(zres(ia(jj)),zu(ia(jj)),zrhs(ia(jj)),
     %      za(ia(jj)),zb(ia(jj)),zc(ia(jj)),nfx,nfy,nfz)

c***
c     we copy all residuals in rhs in order project them to the 
c     next level above
c***

            call copy(zrhs(ia(jj)),zres(ia(jj)),nfx,nfy,nfz)

            nfx=nfx/2+1
            nfy=nfy/2+1
            nfz=nfz/2+1

            call rstrct(zrhs(ia(jj-1)),zrhs(ia(jj)),nfx,nfy,nfz)

c     Initial guess u=0 for the error at the next level above

            call fill0(zu(ia(jj-1)),nfx,nfy,nfz)
      enddo

c      write(*,*) "toto est en bas", NPRE

      call slvsml(zu(ia(1)),zrhs(ia(1)),za(ia(1)),
     %     zb(ia(1)),zc(ia(1)),zei(ia(1)),nnx,nny,nnz,
     %     NPRE,OMEGA)

      do jj=2,NG
         nfx=2*nfx-2
         nfy=2*nfy-2
         nfz=2*nfz-2

c***
c     Use res for temporary storage inside addint, m for storage 
c     of the defect, rhs contains 
c     the previous residuals (identical to 2D version in VMG_relax) 
c 
c     Add interpolation of error u from previous level to this level
c***
         call addint(zu(ia(jj)),zu(ia(jj-1)),zres(ia(jj)),
     %        zm(ia(jj)),zrhs(ia(jj)),za(ia(jj)),zb(ia(jj)),zc(ia(jj)),
     %        nfx,nfy,nfz)

c***
c     relax NPOST times with initial guess e=0
c***
            call fill0(zm(ia(jj)),nfx,nfy,nfz)

            NPOST=RELAX_TIMES*(NG-jj+1)

            do jpost=1,NPOST
               call relax(zm(ia(jj)),zrhs(ia(jj)),za(ia(jj)),
     %              zb(ia(jj)),zc(ia(jj)),zei(ia(jj)),
     %              nfx,nfy,nfz,OMEGA)
            enddo

c     Add correction m  to error u

            call update(zu(ia(jj)),zm(ia(jj)),nfx,nfy,nfz) 

c***
c     calculate residual and copy its value in rhs for next V-cycle.
c***
            call resid(zres(ia(jj)),zm(ia(jj)),zrhs(ia(jj)),
     %      za(ia(jj)),zb(ia(jj)),zc(ia(jj)),nfx,nfy,nfz)

            call copy(zrhs(ia(jj)),zres(ia(jj)),nfx,nfy,nfz)

         enddo

c***        
c     update solution on finest grid id = id + iu
c***
         call update(zd(ia(NG)),zu(ia(NG)),nfx,nfy,nfz)
c         write(*,*) "toto est en haut de nouveau"
         
         norm_old = norm
         norm=h**2*FINDN2S(zrhs(ia(NG)),nfx,nfy,nfz)
         if (printres .eq. 1) then
            write(*,*) "mglin: cycles = ",n,"  res(mean)=",norm
         endif

c***
c     end of V-cycle loop with too many iterations
c***

      enddo
      conv_rate=norm/NORM_OLD
      NORM_OLD=norm
      if (printres .eq. 1) then
         write(*,*) "mglin: return with max cycles="
      endif
  
c***
c     update pressure field
c***

      call copy(zu(ia(NG)),zd(ia(NG)),nfx,nfy,nfz)
      call bc_scalar(zu(ia(NG)),nfx,nfy,nfz)

      cycles = n
      return
      end
c***********************************************************************
      SUBROUTINE rstrct(uc,uf,ncx,ncy,ncz)
c***
      INCLUDE 'undefined.h'
      INTEGER ncx,ncy,ncz,ic,if,jc,jf,kc,kf
      DOUBLE PRECISION uc(ncx,ncy,ncz),uf(2*ncx-2,2*ncy-2,2*ncz-2)
c     SUBROUTINES: bc_scalar
c***
c     each point inside the cube equal 1/8 of the sum of the vertices
c***
      do kc=2,ncz-1
         kf=2*kc-2
         do jc=2,ncy-1
            jf=2*jc-2
            do ic=2,ncx-1
               if = 2*ic - 2
               uc(ic,jc,kc)= .125d0*( uf(if,jf,kf) + uf(if+1,jf,kf) 
     %            + uf(if,jf+1,kf) + uf(if+1,jf+1,kf)
     %            + uf(if,jf,kf+1) + uf(if+1,jf,kf+1)
     %            + uf(if,jf+1,kf+1) + uf(if+1,jf+1,kf+1) )
            enddo
         enddo
      enddo
c***
      call bc_scalar(uc,ncx,ncy,ncz)
c***
      return
      end
c***********************************************************************
      SUBROUTINE interp(uf,uc,nfx,nfy,nfz)
c***
      INCLUDE 'undefined.h'
      INTEGER nfx,nfy,nfz
      DOUBLE PRECISION uc(nfx/2+1,nfy/2+1,nfz/2+1),uf(nfx,nfy,nfz)
      INTEGER ic,if,jc,jf,kc,kf,ncx,ncy,ncz
c     SUBROUTINES: bc_scalar
c***
c     each point of the cube equal to the point inside
c***
      ncx=nfx/2 + 1
      ncy=nfy/2 + 1
      ncz=nfz/2 + 1
      do kc=2,ncz-1
         kf=2*kc-2
         do jc=2,ncy-1
            jf=2*jc-2
            do ic=2,ncx-1
               if=2*ic-2
               uf(if,jf,kf) = uc(ic,jc,kc)
               uf(if,jf+1,kf) = uc(ic,jc,kc)
               uf(if+1,jf,kf) = uc(ic,jc,kc)
               uf(if+1,jf+1,kf) = uc(ic,jc,kc)
               uf(if,jf,kf+1) = uc(ic,jc,kc)
               uf(if,jf+1,kf+1) = uc(ic,jc,kc)
               uf(if+1,jf,kf+1) = uc(ic,jc,kc)
               uf(if+1,jf+1,kf+1) = uc(ic,jc,kc)
            enddo
         enddo
      enddo
c***
      call bc_scalar(uf,nfx,nfy,nfz)
c***
      return
      end
c***********************************************************************
      SUBROUTINE addint(uf,uc,res,error,rhs,a,b,c,nfx,nfy,nfz)
c***
      INCLUDE 'undefined.h'
      INTEGER nfx,nfy,nfz,i,j,k
      DOUBLE PRECISION res(nfx,nfy,nfz),uc(nfx/2+1,nfy/2+1,nfz/2+1),
     %     uf(nfx,nfy,nfz),
     %     a(nfx,nfy,nfz),b(nfx,nfy,nfz),c(nfx,nfy,nfz),
     %     error(nfx,nfy,nfz),rhs(nfx,nfy,nfz)
c     SUBROUTINES: interp
c***
c     interpolate the error from coarse to fine grid
c***
      call interp(error,uc,nfx,nfy,nfz)
c***
c     calculate the new residual with res(new)=res(old)-A*error
c***
      call resid(res,error,rhs,a,b,c,nfx,nfy,nfz)
      call copy(rhs,res,nfx,nfy,nfz)
c***
c correct e(ngrid) with error = Interp(e(ngrid-1))
c***
      do k=1,nfz
         do j=1,nfy
            do i=1,nfx
               uf(i,j,k)=uf(i,j,k)+error(i,j,k)
            enddo
         enddo
      enddo
c***
      return
      end
c***********************************************************************
      SUBROUTINE slvsml(u,rhs,a,b,c,ei,nx,ny,nz,NTOP,OMEGA)
c***
      include 'undefined.h'
      INTEGER NTOP,j,nx,ny,nz
      DOUBLE PRECISION a(nx,ny,nz),b(nx,ny,nz),c(nx,ny,nz)
      DOUBLE PRECISION rhs(nx,ny,nz), u(nx,ny,nz), 
     %     OMEGA,ei(nx,ny,nz)
c     SUBROUTINES: relax
c***
c     call NTOP times relax at the lowest grid
c***  
      do j=1,NTOP
         call relax(u,rhs,a,b,c,ei,nx,ny,nz,OMEGA)
      enddo
      
c***
      return
      end
c***********************************************************************
      SUBROUTINE relax(u,rhs,a,b,c,ei,nx,ny,nz,OMEGA)
c***
      INCLUDE 'undefined.h'
      INTEGER nx,ny,nz, i, ipass, isw, j, jsw, k, ksw
      DOUBLE PRECISION rhs(nx,ny,nz),u(nx,ny,nz)
      DOUBLE PRECISION a(nx,ny,nz),b(nx,ny,nz),c(nx,ny,nz),
     %     ei(nx,ny,nz)
      DOUBLE PRECISION h1, h2, OMEGA, temp
c     SUBROUTINES: bc_scalar
c***

c     estimate the total number of operations. 

c      ZOP = ZOP + nx*ny*nz

      h1=1.d0/(nx-2)
      h2=h1*h1
      ksw=1
c***
      do ipass=1,2
c***
      call bc_scalar(u,nx,ny,nz)
c***
c     solve pressure equation with an alternate sweep
c***
         jsw=ksw
         do k=2,nz-1
           isw=jsw
           do j=2,ny-1
              do i=isw+1,nx-1,2
                 temp = (a(i,j,k)*u(i+1,j,k)+ a(i-1,j,k)*u(i-1,j,k)
     %                  + b(i,j,k)*u(i,j+1,k) + b(i,j-1,k)*u(i,j-1,k)
     %                  + c(i,j,k)*u(i,j,k+1) + c(i,j,k-1)*u(i,j,k-1)
     %                  - h2*rhs(i,j,k))*ei(i,j,k)
                 u(i,j,k) = OMEGA*temp + (1.d0 - OMEGA)*u(i,j,k)
              enddo
              isw=3-isw
           enddo
           jsw=3-jsw
         enddo
         ksw=3-ksw
      enddo
c***
      call bc_scalar(u,nx,ny,nz)
c***
      return
      end
c***********************************************************************
      SUBROUTINE resid(res,u,rhs,a,b,c,nx,ny,nz)
c***
      INCLUDE 'undefined.h'
      INTEGER nx,ny,nz,i,j,k
      DOUBLE PRECISION res(nx,ny,nz),rhs(nx,ny,nz),u(nx,ny,nz)
      DOUBLE PRECISION a(nx,ny,nz),b(nx,ny,nz),c(nx,ny,nz)
      DOUBLE PRECISION h1,h2i
c***
c     calculate the residue of the pressure equation
c***
      h1=1.d0/(nx-2)
      h2i=1.d0/(h1*h1)
      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               res(i,j,k)= - h2i*( a(i,j,k)*(u(i+1,j,k)-u(i,j,k)) 
     %                     + a(i-1,j,k)*(u(i-1,j,k)-u(i,j,k))    
     %                     + b(i,j,k)*(u(i,j+1,k)-u(i,j,k))     
     %                     + b(i,j-1,k)*(u(i,j-1,k) - u(i,j,k))
     %                     + c(i,j,k)*(u(i,j,k+1)-u(i,j,k))   
     %                     + c(i,j,k-1)*(u(i,j,k-1) - u(i,j,k))
     %                     ) + rhs(i,j,k)                    
            enddo
         enddo
      enddo
c***
      call bc_resid(res,nx,ny,nz)
c***
      return
      end
c***********************************************************************
      SUBROUTINE acfromac(a,b,c,ei,a2,b2,c2,nx,ny,nz)
c***
      INCLUDE 'undefined.h'
      INTEGER nx,ny,nz,i,j,k
      DOUBLE PRECISION a(nx,ny,nz),b(nx,ny,nz)
      DOUBLE PRECISION c(nx,ny,nz),ei(nx,ny,nz)
      DOUBLE PRECISION a2(2*nx-2,2*ny-2,2*nz-2)
      DOUBLE PRECISION b2(2*nx-2,2*ny-2,2*nz-2)
      DOUBLE PRECISION c2(2*nx-2,2*ny-2,2*nz-2)
c***
c     with this numerics, a is on the E face, b on the N face, 
c     c on the top, ei in the center of the cell
c***
      do k=2,nz-1
         do j=2,ny-1
            do i=1,nx-1
               a(i,j,k) = 0.25d0*(a2(2*i-1,2*j-1,2*k-1) + 
     %              a2(2*i-1,2*j-2,2*k-1) + a2(2*i-1,2*j-1,2*k-2) +
     %              a2(2*i-1,2*j-2,2*k-2))
            enddo
         enddo
      enddo
      do k=2,nz-1
         do j=1,ny-1
            do i=2,nx-1
               b(i,j,k) = 0.25d0*(b2(2*i-1,2*j-1,2*k-1) +
     %              b2(2*i-2,2*j-1,2*k-1) + b2(2*i-1,2*j-1,2*k-2) +
     %              b2(2*i-2,2*j-1,2*k-2))
            enddo
         enddo
      enddo
      do k=1,nz-1
         do j=2,ny-1
            do i=2,nx-1
               c(i,j,k) = 0.25d0*(c2(2*i-1,2*j-1,2*k-1) +
     %              c2(2*i-2,2*j-1,2*k-1) + c2(2*i-1,2*j-2,2*k-1) +
     %              c2(2*i-2,2*j-2,2*k-1))
            enddo
         enddo
      enddo
c***
      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               ei(i,j,k) = 1.d0/(a(i,j,k) + a(i-1,j,k) + b(i,j,k) +
     %                           b(i,j-1,k) + c(i,j,k) + c(i,j,k-1))
            enddo
         enddo
      enddo
c***
      return
      end
c***********************************************************************
      SUBROUTINE update(aout,ain,nx,ny,nz)
c***
      include 'undefined.h'
      INTEGER nx,ny,nz
      DOUBLE PRECISION ain(nx,ny,nz),aout(nx,ny,nz)
      
      INTEGER i,j,k
      do k=2,nz-1
         do j =2,ny-1
            do i=2,nx-1
               aout(i,j,k)=aout(i,j,k)+ain(i,j,k)
            enddo
         enddo
      enddo

c***
      call bc_scalar(aout,nx,ny,nz)
c***
      return
      end