
      SUBROUTINE filter(c, cfilt, nx,ny)
      include 'undefined.h'
      INTEGER nx,ny,i,j
      DOUBLE PRECISION c(nx,ny), cfilt(nx,ny) 
      do j=2,ny-1
         do i=2,nx-1
            cfilt(i,j) = 0.25*c(i,j) 
     *  + 0.125*(c(i-1,j) + c(i+1,j) + c(i,j-1) + c(i,j+1))
     *  + 0.0625*(c(i-1,j-1) + c(i+1,j-1) + c(i-1,j+1) + c(i+1,j+1))
         enddo
      enddo
      return
      end

      DOUBLE PRECISION FUNCTION normmean(u,nx,ny)
      include 'undefined.h'
      integer  nx,ny, i, j 
      DOUBLE PRECISION U(NX,NY)
      normmean=0.0d0
      do j=2,ny-1
         do i=2,nx-1
            normmean = normmean + (u(i,j))**2
         enddo
      enddo
      normmean = normmean*(1.0d0/(nx-2))*(1.0d0/(ny-2))
      normmean = dsqrt(normmean)
      return
      end

      DOUBLE PRECISION FUNCTION xo(u,nx,ny)
            include 'undefined.h'
      integer  nx,ny, i, j 
      DOUBLE PRECISION U(NX,NY)
      xo=u(1,2)
      return
      end
      DOUBLE PRECISION FUNCTION yo(u,nx,ny)
            include 'undefined.h'
      integer  nx,ny, i, j 
      DOUBLE PRECISION U(NX,NY)
      yo=u(2,1)
      return
      end

      
      DOUBLE PRECISION FUNCTION findmax(u,nx,ny,imax,jmax)
      include 'undefined.h'
      integer imax, jmax, nx,ny, i, j 
      DOUBLE PRECISION U(NX,NY),  res
      findmax=-1.d10
      do j=2,ny-1
         do i=2,nx-1
            res = u(i,j)
            if(res.gt.findmax) then
               findmax = res
               imax = i
               jmax = j
            endif
         enddo
      enddo
      return
      end

      DOUBLE PRECISION FUNCTION momc(v,cc,r1,nx,ny)
      include 'undefined.h'
      integer i,j,nx,ny
      DOUBLE PRECISION v(nx,ny),cc(nx,ny),r1,r2

      momc=0.d0

      r2=1.d0
      do j=2,ny-1
        do i=2,nx-1
      momc=momc+
     * 0.5* (i-1.5d0)*(r2*(cc(i,j)) + r1*(1-cc(i,j)))*(v(i,j)+v(i,j+1))

        enddo
      enddo

      return
      end

      DOUBLE PRECISION FUNCTION momc2(v,cc,r1,nx,ny)
      include 'undefined.h'
      integer i,j,nx,ny
      DOUBLE PRECISION v(nx,ny),cc(nx,ny),r1,r2,vol

      momc2=0.d0
      vol=0.d0

      r2=(1.d0-r1)

      do j=2,ny-1
        do i=2,nx-1
      momc2=momc2+
     * 0.5*(i-1.5d0)*(cc(i,j))*(v(i,j)+v(i,j+1))
      vol=vol + (i-1.5d0)*cc(i,j)
        enddo
      enddo
      momc2=momc2/vol
      return
      end

      DOUBLE PRECISION FUNCTION fmodmax1(u,nx,ny,imax,jmax)
      include 'undefined.h'
      integer imax, jmax, nx,ny, i, j 
      DOUBLE PRECISION U(NX,NY),  res
      fmodmax1=0.d0
      do j=3,ny-1
         do i=3,nx-1
c            res = u(i,j)/(i-1.5d0)
            res = u(i,j)
            if(dabs(res).gt.fmodmax1) then
               fmodmax1 = dabs(res)
               imax = i
               jmax = j
            endif
         enddo
      enddo
      return
      end
      DOUBLE PRECISION FUNCTION fmodmax(u,nx,ny,imax,jmax)
      include 'undefined.h'
      integer imax, jmax, nx,ny, i, j 
      DOUBLE PRECISION U(NX,NY),  res
      fmodmax=0.d0
      do j=3,ny-1
         do i=3,nx-1
            res = u(i,j)
            if(dabs(res).gt.fmodmax) then
               fmodmax = dabs(res)
               imax = i
               jmax = j
            endif
         enddo
      enddo
      return
      end

      DOUBLE PRECISION FUNCTION findmin(u,nx,ny,imin,jmin)
      include 'undefined.h'
      integer imin, jmin, nx,ny, i, j 
      DOUBLE PRECISION U(NX,NY),  res
      findmin=1.d100
      do j=2,ny-1
         do i=2,nx-1
            res = u(i,j)
            if(res.lt.findmin) then
               findmin = res
               imin = i
               jmin = j
            endif
         enddo
      enddo
      return
      end

      SUBROUTINE finddiff(u,v,y,nx,ny,imax,jmax)
      include 'undefined.h'
      integer imax, jmax, nx,ny, i, j 
      DOUBLE PRECISION U(NX,NY), v(nx,ny) , y, res
      y=0.d0
      
      do j=2,ny-1
         do i=2,nx-1
            res = u(i,j)  - v(i,j)
            if(dabs(res).gt.y) then
               y = dabs(res)
               imax = i
               jmax = j
            endif
         enddo
      enddo
      return
      end

      DOUBLE PRECISION FUNCTION sumfield(u,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny,i,j
      DOUBLE PRECISION u(nx,ny)
c
      sumfield=0.d0
      do j=2,ny-1
         do i=2,nx-1
c            write(6,*) u(i,j),i,j
            sumfield = sumfield + u(i,j) 
         enddo
      enddo

      return
      end

      DOUBLE PRECISION FUNCTION sumfield2(u,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny,i,j
      DOUBLE PRECISION u(nx,ny)

      sumfield2=0.d0
      do j=2,ny-1
         do i=2,nx-1
c            write(53,*) u(i,j),i,j
            sumfield2 = sumfield2 + (i-1.5d0)*u(i,j) 
         enddo
      enddo
      return
      end

      DOUBLE PRECISION FUNCTION sumfield3(nx,ny)
      include 'undefined.h'
      INTEGER nx,ny,i,j

      sumfield3=0.d0
      do j=2,ny-1
         do i=2,nx-1
            sumfield3 = sumfield3 + (i-1.5d0) 
         enddo
      enddo
      return
      end

      DOUBLE PRECISION FUNCTION usang(tau,v,c,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny,i,j
      DOUBLE PRECISION v(nx,ny),c(nx,ny),voltot,volliq,volout
      DOUBLE PRECISION uliq,uout,tau,h
      h=1.0d0/(nx-2)
      usang=0.0d0
      volliq=0.0d0
      volout=0.0d0
      voltot=0.d0
      uliq=0.0d0
      uout=0.0d0
      do j=2,ny-1
         do i=2,nx-1
            volliq = volliq + c(i,j)*(i-1.5d0) 
            voltot = voltot + (i-1.5d0) 
            uliq = uliq + c(i,j)*0.5d0*(i-1.5d0)*(v(i,j)+v(i,j+1))
            uout = uout +(1.d0-c(i,j))*0.5d0*(i-1.5d0)*(v(i,j)+v(i,j+1))
         enddo
      enddo
      uliq=uliq*(1.d0/volliq)
      volout=voltot-volliq
      uout=uout*(1.d0/volout)
      usang=(h/tau)*(volout/voltot)*(uliq-uout)
      return
      end

      DOUBLE PRECISION FUNCTION findn2(u,v,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny,i,j
      DOUBLE PRECISION u(nx,ny), v(nx,ny)
      findn2=0.d0
      do j=3,ny-1
         do i=2,nx-1
            findn2 = findn2 + u(i,j)**2 + v(i,j)**2
         enddo
      enddo
      findn2 = dsqrt(findn2/(nx*ny))
      return
      END




      DOUBLE PRECISION FUNCTION findninf(u,v,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny,i,j
      DOUBLE PRECISION u(nx,ny), v(nx,ny),vv
      findninf=0.d0
      do j=3,ny-1
         do i=2,nx-1
            vv = dsqrt(u(i,j)**2 + v(i,j)**2)
            findninf=dmax1(vv,findninf)
         enddo
      enddo
      return
      END

      DOUBLE PRECISION FUNCTION findninf1(u,v,c,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny,i,j
      DOUBLE PRECISION u(nx,ny), v(nx,ny),vv,c(nx,ny)
      findninf1=0.d0
      do j=3,ny-1
         do i=2,nx-1
            if (c(i,j) .ge. 0.5d0) then
            vv = dsqrt(u(i,j)**2 + v(i,j)**2)
            findninf1=dmax1(vv,findninf1)
            endif
         enddo
      enddo
      return
      END

      subroutine iniampli(nx,ny)
      include 'undefined.h'
      integer nx,ny
      
      open(unit=17,file="hx.par",form="formatted")
      close(unit=17,status="delete")
      open(unit=17,file="hx.par",form="formatted",status="new")
      open(unit=18,file="hy.par",form="formatted")
      close(unit=18,status="delete")
      open(unit=18,file="hy.par",form="formatted",status="new")
      return
      end

      DOUBLE PRECISION FUNCTION hofy(u,v,nx,ny)
      include 'undefined.h'
      integer nx,ny,record,i,j
      DOUBLE PRECISION u(nx,ny),v(nx,ny)
      common /pvel/record
      include 'machdep.f'
      include 'mach.h'
      
      do j=2,ny-1
         v(1,j)=0.0d0
      enddo
      
      do j=2,ny-1
         do i=2,nx-1
         v(1,j)=v(1,j)+u(i,j)
         enddo
      enddo

      do j=2,ny-1
         write(17,*) v(1,j)
      enddo

      call flush(17)
      return 
      end

      DOUBLE PRECISION FUNCTION hofx(u,v,nx,ny)
      include 'undefined.h'
      integer nx,ny,record,i,j
      DOUBLE PRECISION u(nx,ny),v(nx,ny)
      common /pvel/record
      include 'machdep.f'
      include 'mach.h'
      
      do i=2,nx-1
         v(i,1)=0.0d0
      enddo
      
      do i=2,nx-1
         do j=2,ny-1
            v(i,1)=v(i,1)+u(i,j)
         enddo
      enddo

      do i=2,nx-1
         write(18,*) v(i,1)
      enddo

      call flush(18)
      return 
      end

      DOUBLE PRECISION FUNCTION findmax3(u,nx,ny,imax)
      include 'undefined.h'
      integer imax, nx,ny, i, j 
      DOUBLE PRECISION U(NX,NY),  res
      findmax3=-1.d10
      do i=2,nx-1
         res = u(i,1)
         if(res.gt.findmax3) then
            findmax3 = res
            imax = i
         endif
      enddo
      return
      end

      DOUBLE PRECISION FUNCTION findmax2(u,nx,ny,jmax)
      include 'undefined.h'
      integer jmax, nx,ny, i, j 
      DOUBLE PRECISION U(NX,NY),  res
      findmax2=-1.d10
      do j=2,ny-1
         res = u(1,j)
         if(res.gt.findmax2) then
            findmax2 = res
            jmax = j
         endif
      enddo
      return
      end



