c************************************************************
      SUBROUTINE initpressure(nx,ny)
      include 'undefined.h'
      INTEGER mymalloc,ngx,ngy,jx,jy,ncycle,cycles,npre,npost
      INTEGER j,nn,nx,ny,nnx,nny,ifetch,dfetch,nres,relax_times,
     * relax_times_old,ncycle_old,printres
      DOUBLE PRECISION mineps,maxeps,norm_old,GP_dfetch
      INCLUDE 'size.h'
      INCLUDE 'pointers.f'
      COMMON /cycling/ mineps,maxeps,cycles,ncycle,npre,npost,nres,
     * relax_times,norm_old,printres,relax_times_old,ncycle_old

C     Initialize memory allocation scheme in mymalloc (mymalloc.f) 
      mem=0

      if(ifetch("MG_cycles_bound",cycles).ne.1) 
     *  stop 'missing MG_cycles_bound'
      if(ifetch("nres",nres).ne.1) nres=0 
      if (ifetch("printres",printres).ne.1) printres=0
      j = dfetch("div_tolerance",mineps)
      if(ifetch("relaxfine",relax_times).ne.1) relax_times=3
      ncycle=cycles
      write(6,*) "MG_cycles_bound=",ncycle  
      write(6,*) "div_tolerance=",mineps
      write(6,*) "printres=",printres

c      OMEGA=1.0d0
      OMEGA = GP_dfetch("OMEGA",1.d0,"WARNING")
   
      write(6,*) "OMEGA=",OMEGA
      nn=min(nx,ny)
c
C     compute ng 
c
      j = nn - 2
      jx= nx-2
      jy= ny-2
      ng=1
      ngx=1
      ngy=1
 1    if(j.gt.2) then
         j = j/2
         ng = ng+ 1
         go to 1
      endif
      if((2**ng + 2).ne.nn ) then
         write(6,*) "nn= ",nn,"   ng= ",ng
         stop "initpressure: nn.ne. 2**ng + 2"
      endif

 2    if(jx.ge.4) then 
         jx = jx/2
         ngx = ngx+ 1
         go to 2
      endif
c      if((2**ngx + 2).ne.nx ) then 
c         write(6,*) "nx=",nx,"   ngx=",ngx
c         stop "initpressure: nnx.ne. 2**ngx + 2"
c      endif
 3    if(jy.ge.4) then
         jy = jy/2
         ngy = ngy+ 1
         go to 3
      endif
c      if((2**ngy + 2).ne.ny ) then
c         write(6,*) "ny=",ny,"   ngy=",ngy
c         stop "initpressure: ny.ne. 2**ngy + 2"
c      endif
      nnx=nx
      nny=ny
         iccol = mymalloc(nnx*nny)
      do j=ng,1,-1
         ira(j) = mymalloc(nnx*nny)
         irc(j) = mymalloc(nnx*nny)
         irei(j) = mymalloc(nnx*nny)
         icc(j) = mymalloc(nnx*nny)
         iu(j) = mymalloc(nnx*nny)
         ia(j) = mymalloc(nnx*nny)
         ic(j) = mymalloc(nnx*nny)
         iei(j) = mymalloc(nnx*nny)
         ib(j) = mymalloc(nnx*nny)
         id(j) = mymalloc(nnx*nny)
         ires(j) = mymalloc(nnx*nny)
         irhs(j) = mymalloc(nnx*nny)
         nnx=nnx/2+1
         nny=nny/2+1
         write(6,*) "level ", j, " nx = ",nnx -2 ," ny = ",nny - 2
      enddo
      print*,'MEMLEN,mem:',MEMLEN,mem
      PRINT*,'Gbytes required: ', (MEMLEN*64.d0)/1.d9, (mem*64.d0)/1.d9

      return
      end

c************************************************************
      integer function mymalloc(len)
      include 'undefined.h'
      INTEGER len
      INCLUDE 'size.h'
      INCLUDE 'pointers.f'
      if(mem+len.gt.MEMLEN)  then 
         write(6,*) 'MEMLEN=',MEMLEN,'mem+len',mem+len
         STOP 'insufficient memory in mymalloc'
      endif
      mymalloc = mem +1
      mem = mem + len
      return
      END
c
c
c
c************************************************************
c   This subroutine inverts the pressure equation
c
c   ON INPUT:
c  
c   nx by ny double arrays:
c   ---------------------
c   irhs(ng):  initial right hand side (divergence field)
c   cc: concentrations or volume fractions of phase 1
c   a:  coefficients for equation on grid n, computed with "afromc"
c   c:  coefficients "   "
c   ei: coefficients inverse "  "
c
c   scalar parameters:
c   ------------------
c   r1, r2: parameters used by afromc.
c   nx,ny:      dimension of arrays
c   ncycle: number of cycles in Vcycle
c
c  ON OUTPUT:
c
c   everything the same, except
c
c   iu(ng):      contains the pressure 
c

      SUBROUTINE mglin(r1,r2,nx,ny)
C     uses addint copy fill0 interp relax resid rstrct slvsml afromc
      include 'undefined.h'
      INTEGER nx,ny,nn,nnx,nny,ncycle,cycles,NPRE,NPOST,n,nres,ifetch
      INTEGER i,jj,jpost,jpre,nfx,nfy,ngrid,relax_times,
     *  printres,relax_times_old,ncycle_old
      DOUBLE PRECISION r1, r2
      DOUBLE PRECISION mineps,maxeps,normmean,norm,norm_old,conv_rate
      SAVE conv_rate
      COMMON /cycling/ mineps,maxeps,cycles,ncycle,npre,npost,nres,
     *relax_times,norm_old,printres,relax_times_old,ncycle_old
      INCLUDE 'size.h'
      INCLUDE 'pointers.f'
      DATA conv_rate /0.5d0/
      nn=min(nx,ny)
      if((2**ng+2).ne.nn) then
           write(6,*) "nn=",nn,"nx=",nx,"ny=",ny,"ng=",ng
           stop "fmg: nn.ne.(2**ng+2)"
      endif

      i=ifetch("MG_cycles_bound",ncycle)
      ngrid=NG
      nnx=nx
      nny=ny
      nfx=nnx
      nfy=nny
      call coli(z(ira(ngrid)),z(irc(ngrid)),z(irei(ngrid)),z(ia(ngrid))
     *   ,z(ic(ngrid)),nnx,nny)
 1    if (nn.gt.4) then
         nn=nn/2 + 1
         nnx=nnx/2 + 1
         nny=nny/2+1
         ngrid = ngrid - 1
         
         call acfromac(z(ira(ngrid)),z(irc(ngrid)),
     *        z(irei(ngrid)),z(ira(ngrid+1)),z(irc(ngrid+1)),nnx,nny)
         go to 1
      endif
      
      
c     first copy previous solution at time t^n in a field d(NG)
      call copy(z(id(NG)),z(iu(NG)),nfx,nfy)
      
c     then 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.
      
      call resid(z(ires(NG)),z(id(NG)),z(irhs(NG)),z(ira(NG)),
     *     z(irc(NG)),nfx,nfy)
      call copy(z(irhs(NG)),z(ires(NG)),nfx,nfy)
      norm=normmean(z(irhs(NG)),nfx,nfy)
      if (printres .eq. 1) then 
         write(6,*) "beginning MG"," res(mean)=",
     *        norm," max number of cycles=",ncycle
      endif
c     if bad convergence rate at the end of last MG cycle, then increase the number of smoothing 
c     for the new cycle
      relax_times_old=relax_times
      ncycle_old=cycles
      if (conv_rate .gt. 0.6) then 
         relax_times=relax_times + 1
      else if (conv_rate .lt. 0.1)then
         relax_times=max(relax_times-1,1)
      endif
      if (printres .eq. 1) then 
      write(6,*) "numb relax on fine grid=",relax_times," conv_rate ="
     *        ,conv_rate,"norm_old=",norm_old
      endif
      do n=1,ncycle
c     write(6,*) "numero cycle",n
         if(norm.lt.mineps.and.n.gt.1) then
            call copy(z(iu(NG)),z(id(NG)),nfx,nfy)
            call makebcpf(z(iu(NG)),nfx,nfy)
            cycles=n-1
            if (printres .eq. 1) then 
               write(6,*) "mglin: return cycles=",cycles,"norm=",norm
            endif
            conv_rate=norm/norm_old
            return
         endif
      if(printres.eq.1) then
         write(6,*) "debut V-cycle","   res(mean)=",
     * normmean(z(irhs(NG)),nfx,nfy)
         write(6,*) "norm before old", norm
      endif 
      norm_old=norm

c     on fine grid relax on the error e(NG) with initial guess e(NG)=0 
         call fill0(z(iu(NG)),nfx,nfy)
c     Now u contains the error (or correction) and not the solution anymore.
         do jj=NG,2,-1
c     compute 1 relaxation on the NGth grid, 2 on the NG-1th ... NG-1 on
c     the 2nd (if relax_times=1)
            NPRE=relax_times*(NG-jj+1)
            do jpre=1,NPRE
               call relax(z(iu(jj)),z(irhs(jj)),
     *              z(ira(jj)),z(irc(jj)),z(irei(jj)),nfx,nfy,OMEGA)
c     we copy all residuals in rhs in order to use it later for the correction step 
c     (in addint)
            enddo
            call resid(z(ires(jj)),z(iu(jj)),z(irhs(jj)),
     *           z(ira(jj)),z(irc(jj)),nfx,nfy)
            call copy(z(irhs(jj)),z(ires(jj)),nfx,nfy)
c     write(6,*) "descente relax num.",jpre,"res(mean)=",
c     * normmean(z(ires(jj)),nfx,nfy)
                        
            nfx=nfx/2+1
            nfy=nfy/2+1
c     write(6,*) "toto descend"
            call rstrct(z(irhs(jj-1)),z(irhs(jj)),nfx,nfy)
            call fill0(z(iu(jj-1)),nfx,nfy)
         enddo
c     write(6,*) "toto est tout en bas"
         call slvsml(z(iu(1)),z(irhs(1)),
     *        z(ira(1)),z(irc(1)),z(irei(1)),nfx,nfy,OMEGA)
         
         do jj=2,NG
            nfx=2*nfx-2
            nfy=2*nfy-2
c     Use res and temporary storage inside addint, b for storage of the defect, rhs contains 
c     the previous residuals
            call addint(z(iu(jj)),z(iu(jj-1)),z(ires(jj)),z(ib(jj)),
     *           z(irhs(jj)),z(ira(jj)),z(irc(jj)),nfx,nfy)


c relax NPOST times with initial guess e=0
            call fill0(z(ib(jj)),nfx,nfy)
            NPOST=relax_times*(NG-jj+1)
            do jpost=1,NPOST
               call relax(z(ib(jj)),z(irhs(jj)),
     *              z(ira(jj)),z(irc(jj)),z(irei(jj)),nfx,nfy,OMEGA)
c     calculate residual and copy its value in rhs for next V-cycle.
            enddo
            call update(z(iu(jj)),z(ib(jj)),nfx,nfy) 
            call resid(z(ires(jj)),z(ib(jj)),z(irhs(jj)),
     *           z(ira(jj)),z(irc(jj)),nfx,nfy)
            call copy(z(irhs(jj)),z(ires(jj)),nfx,nfy)
c     write(6,*) "montee relax num.",jpost,"res(mean)=",
c     *           normmean(z(iu(jj)),nfx,nfy)
            
         enddo
         
c     update solution on finest grid id = id+iu
         call update(z(id(NG)),z(iu(NG)),nfx,nfy)
c     write(6,*) "toto est en haut de nouveau"
         norm=normmean(z(irhs(NG)),nfx,nfy)
         if (printres .eq. 1) then 
            if (nres .eq. 1) then
               write(6,*) " ncycle=",n,"res(mean)=",norm
            endif
         endif
         conv_rate=norm/norm_old
         if(printres.eq.1) then 
            write(6,*) "norm before conv-rate",norm
            write(6,*) "conv_rate", conv_rate
         endif
c end of V-cycle loop
      enddo
      if (printres .eq. 1) then 
         write(6,*) "mglin: MG_cycles_bound is reached, norm=",norm
      endif
c if bad convergence rate : increase the number of relaxations
         if (conv_rate .gt. 0.6) then 
            relax_times=relax_times + 1
         else if (conv_rate .lt. 0.1) then
            relax_times=max(relax_times-1,1)
         endif
c****** debug
c               write(6,*) "conv_rate", conv_rate
c         write(6,*) "relax_fine_grid=",relax_times,"conv_rate_last_cycle ="
c     *    ,conv_rate
c     update pressure field
      call copy(z(iu(NG)),z(id(NG)),nfx,nfy)
c **DG to verify
      call makebcpf(z(iu(NG)),nfx,nfy)
c **
      return
      end


c
c***********************************************************************
c
      SUBROUTINE rstrct(uc,uf,ncx,ncy)
      include 'undefined.h'
      INTEGER ncx,ncy
      DOUBLE PRECISION uc(ncx,ncy),uf(2*ncx-2,2*ncy-2)
      INTEGER ic,if,jc,jf
      DO jc=2,ncy-1
         jf=2*jc-2
         do ic=2,ncx-1
            if = 2*ic - 2
            uc(ic,jc)= .25d0*(uf(if,jf) + uf(if+1,jf) 
     *         + uf(if,jf+1) + uf(if+1,jf+1))
         enddo
      enddo

      call makebcpf(uc,ncx,ncy)

      return
      END


c
c***********************************************************************
c
      SUBROUTINE interp(uf,uc,nfx,nfy)
      include 'undefined.h'
      INTEGER nfx,nfy
      DOUBLE PRECISION uc(nfx/2+1,nfy/2+1),uf(nfx,nfy)
      INTEGER ic,if,jc,jf,ncx,ncy
      
      ncx=nfx/2 + 1
      ncy=nfy/2 + 1
      do jc=2,ncy-1
         jf=2*jc-2
         do ic=2,ncx-1
            if=2*ic-2
            uf(if,jf) = uc(ic,jc)
            uf(if,jf+1) = uc(ic,jc)
            uf(if+1,jf) = uc(ic,jc)
            uf(if+1,jf+1) = uc(ic,jc)
c            uf(if,jf) =.0625d0*(uc(ic-1,jc-1)+3.0d0*uc(ic-1,jc)
c     &                  +3.0d0*uc(ic,jc-1)+9.0d0*uc(ic,jc))
c            uf(if,jf+1)=.0625d0*(uc(ic-1,jc+1)+3.0d0*uc(ic-1,jc)
c     &                    +3.0d0*uc(ic,jc+1)+9.0d0*uc(ic,jc))
c            uf(if+1,jf)=.0625d0*(uc(ic+1,jc-1)+3.0d0*uc(ic,jc-1)
c     &                   +3.0d0*uc(ic+1,jc)+9.0d0*uc(ic,jc))
c            uf(if+1,jf+1)=.0625d0*(uc(ic+1,jc+1)+3.0d0*uc(ic+1,jc)
c     &                   +3.0d0*uc(ic,jc+1)+9.0d0*uc(ic,jc))
         enddo
      enddo

      call makebcpf(uf,nfx,nfy)
      
      return
      END

c
c***********************************************************************
c
      SUBROUTINE addint(uf,uc,res,error,rhs,a,c,nfx,nfy)
      include 'undefined.h'
      INTEGER nfx,nfy
      DOUBLE PRECISION res(nfx,nfy),uc(nfx/2+1,nfy/2+1),uf(nfx,nfy)
      DOUBLE PRECISION error(nfx,nfy),rhs(nfx,nfy),a(nfx,nfy),c(nfx,nfy)

      INTEGER i,j
c  interpolate the error from coarse to fine grid
      call interp(error,uc,nfx,nfy)
c calculate the new residual with res(new)=res(old)-A*error
      call resid(res,error,rhs,a,c,nfx,nfy)
      call copy(rhs,res,nfx,nfy)
c correct e(ngrid) with error = Interp(e(ngrid-1))
      do j=2,nfy-1
         do i=2,nfx-1
            uf(i,j)=uf(i,j)+error(i,j)
         enddo
      enddo
       call makebcpf(uf,nfx,nfy)
      
      return
      END

c
c***********************************************************************
c
      SUBROUTINE slvsml(u,rhs,a,c,ei,nx,ny,OMEGA)
      include 'undefined.h'
      INTEGER ntop
      INTEGER j,nx,ny
      DOUBLE PRECISION a(nx,ny),c(nx,ny),ei(nx,ny)
      DOUBLE PRECISION rhs(nx,ny),u(nx,ny),OMEGA
c The system is solved "exaclty" on coarsest grid using a great number of 
c relaxations : the work is small since the number of points is small.
      NTOP=50
      do j=1,NTOP
         call relax(u,rhs,a,c,ei,nx,ny,OMEGA)
      enddo
      return
      END
c
c********************************************************************
***
c
      SUBROUTINE relax(u,rhs,a,c,ei,nx,ny,OMEGA)
      include 'undefined.h'
      INTEGER nx,ny
      DOUBLE PRECISION rhs(nx,ny),u(nx,ny),OMEGA,temp
      DOUBLE PRECISION a(nx,ny),c(nx,ny),ei(nx,ny)
      INTEGER i, ipass, isw, j, jsw
      DOUBLE PRECISION h, h2
      h=1.d0/(nx-2)
      h2=h*h
      jsw=1
c
c Force periodicity in x direction
c:
      do ipass=1,2
         call makebcpf(u,nx,ny)
         isw=jsw
         do j=2,ny-1
            do i=isw+1,nx-1,2
               temp = (a(i,j)*u(i+1,j)+ a(i-1,j)*u(i-1,j)
     *                + c(i,j)*u(i,j+1) + c(i,j-1)*u(i,j-1)
     *                - h2*rhs(i,j))*ei(i,j)
               u(i,j) = OMEGA*temp + (1.0d0 - OMEGA)*u(i,j)
            enddo
            isw=3-isw
         enddo
         jsw=3-jsw
      enddo

      call makebcpf(u,nx,ny)


      return
      END

c
c***********************************************************************
c

      SUBROUTINE resid(res,u,rhs,a,c,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny
      DOUBLE PRECISION res(nx,ny),rhs(nx,ny),u(nx,ny)
      DOUBLE PRECISION a(nx,ny),c(nx,ny)
      INTEGER i,j
      DOUBLE PRECISION h,h2i,alpha
      h=1.d0/(nx-2)
      h2i=1.d0/(h*h)
      alpha=0.0d0
      do j=2,ny-1
         do i=2,nx-1
            res(i,j)= - h2i*(a(i,j)*(u(i+1,j)-u(i,j))
     *                +a(i-1,j)*(u(i-1,j)-u(i,j))
     *                + c(i,j)*(u(i,j+1)-u(i,j)) 
     *                + c(i,j-1)*(u(i,j-1) - u(i,j))
     *                ) + rhs(i,j)
c            alpha=alpha+res(i,j)
         enddo
      enddo
c
c
c      write(6,*) "sum res =",alpha
      call makebcresf(res,nx,ny)

      return
      END

c
c***********************************************************************
c
      SUBROUTINE copy(aout,ain,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny
      DOUBLE PRECISION ain(nx,ny),aout(nx,ny)
      
      INTEGER i,j
      do j =1,ny
         do i=1,nx
            aout(i,j)=ain(i,j)
         enddo
      enddo
      return
      END

c
c***********************************************************************
c
      SUBROUTINE fill0(u,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny
      DOUBLE PRECISION u(nx,ny)
      
      INTEGER i,j
      do j=1,ny
         do i=1,nx
            u(i,j)=0.d0
         enddo
      enddo
      return
      END


c
c***********************************************************************
c
      SUBROUTINE store(u,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny
      DOUBLE PRECISION u(nx,ny)

      INTEGER i,j
      print*,"store,nx,ny=",nx,ny
      do j=1,ny
         do i=1,nx
            write (6,*) u(i,j)
         enddo
      enddo
      stop 'end store'
      return
      END






c***********************************************************************
c
      SUBROUTINE acfromac(a,c,ei,a2,c2,nx,ny)
      include 'undefined.h'
      INTEGER i,j,nx,ny
      DOUBLE PRECISION a(nx,ny),c(nx,ny),ei(nx,ny)
      DOUBLE PRECISION a2(2*nx-2,2*ny-2),c2(2*nx-2,2*ny-2)

         do j=2,ny-1
            do i=1,nx-1
               a(i,j) = 0.5d0*(a2(2*i-1,2*j-2)+a2(2*i-1,2*j-1))
            enddo
         enddo
      
         do j=1,ny-1
            do i=2,nx-1
               c(i,j) = 0.5d0*(c2(2*i-2,2*j-1)+c2(2*i-1,2*j-1)) 
            enddo 
         enddo 
         do j=2,ny-1
            do i=2,nx-1
               ei(i,j) = 1.0d0/(a(i-1,j)+a(i,j)+c(i,j-1)+c(i,j))
         enddo 
      enddo 
      return
      END
c
c***********************************************************************
c

      SUBROUTINE eifromac(a,c,ei,r1,r2,nx,ny)
      include 'undefined.h' 
      INTEGER nx,ny 
      DOUBLE PRECISION a(nx,ny),c(nx,ny),ei(nx,ny) 
      DOUBLE PRECISION r1,r2
      INTEGER i,j 
 
      do j=2,ny-1 
         do i=2,nx-1 
             a(i,j) = 1.d0/(r1*a(i,j) + r2) 
             c(i,j) = 1.d0/(r1*c(i,j) + r2) 
          enddo 
      enddo 
 
      do j=2,ny-1 
          a(1,j) = a(2,j) 
      enddo 
 
      do j=2,nx-1 
c          c(j,1) =  1.d0/(r1*c(j,1) + r2) 
          c(j,1) =  c(j,2) 
c          c(j,ny-1) = 0.0d0
      enddo

      do j=2,ny-1
         do i=2,nx-1
             ei(i,j) = 1.d0/(a(i,j) + a(i-1,j) + c(i,j) + c(i,j-1))
         enddo
      enddo

      return
      end
 

c
c***********************************************************************
c
      SUBROUTINE update(aout,ain,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny
      DOUBLE PRECISION ain(nx,ny),aout(nx,ny)
      
      INTEGER i,j
      do j =2,ny-1
         do i=2,nx-1
            aout(i,j)=aout(i,j)+ain(i,j)
         enddo
      enddo
      
      call makebcpf(aout,nx,ny)
      return
      END

c
c***********************************************************************
c
      SUBROUTINE writefield(u,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny
      DOUBLE PRECISION u(nx,ny)
      
      INTEGER i,j
      do j =1,ny
         do i=1,nx
            write(6,*) "p",i,",",j,"=",u(i,j)
         enddo
      enddo
      return
      END


c
c***********************************************************************
c

      SUBROUTINE coli(re,rf,rei,e,f,nx,ny)
      include 'undefined.h'
      INTEGER i,j,nx,ny
      DOUBLE PRECISION e(nx,ny),f(nx,ny),h
      DOUBLE PRECISION re(nx,ny),rf(nx,ny),rei(nx,ny)
      h=1.0d0/(nx-2)
c
c     e,f are oringinal a,c fields
c     re, rf are a, c fields multiplied by radius r (factor h removed)
c
c
c     rf(1,j) is never used
c
      do j=1,ny
         do i=1,nx
            re(i,j)=(i-1)*e(i,j)
            rf(i,j)=(i-1.5d0)*f(i,j)
         enddo
      enddo
c
c     For inlet boundary conditions, annul pressure correction for
c      v(i,j=2) -> annul c(i,1) 

c     do i=1,nx
c        rf(i,2)=0.d0
c      enddo
c
      do j=2,ny
         do i=2,nx
            rei(i,j)=1.0d0/(re(i,j)+re(i-1,j)+rf(i,j)+rf(i,j-1))
        enddo
      enddo

      return
      END
c
c***********************************************************************
c


