c***********************************************************************
      DOUBLE PRECISION FUNCTION FINDMAX(u,nx,ny,nz,imax,jmax,kmax)
c***
c     get the max value, and its indices, of a scalar field
c***
      INCLUDE 'undefined.h'
      INTEGER imax, jmax, kmax, nx, ny, nz, i, j, k 
      DOUBLE PRECISION u(nx,ny,nz)
c***
      FINDMAX = u(2,2,2)
      imax = 2
      jmax = 2
      kmax = 2
      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               if (u(i,j,k) .gt. FINDMAX) then
                  FINDMAX = u(i,j,k)
                  imax = i
                  jmax = j
                  kmax = k
               endif
            enddo
         enddo
      enddo
c***
      return
      end
c***********************************************************************
      DOUBLE PRECISION FUNCTION FMODMAX(u,nx,ny,nz,imax,jmax,kmax)
c***
c     get the absolute max value, and its indices, of a scalar field
c***
      INCLUDE 'undefined.h'
      INTEGER imax, jmax, kmax, nx, ny, nz, i, j, k 
      DOUBLE PRECISION u(nx,ny,nz), res
      INTRINSIC DABS
c***
      FMODMAX=-1.d0
      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               res = DABS(u(i,j,k))
               if (res .gt. FMODMAX) then
                  FMODMAX = res
                  imax = i
                  jmax = j
                  kmax = k
               endif
            enddo
         enddo
      enddo
c***
      return
      end
c***********************************************************************
      DOUBLE PRECISION FUNCTION FINDMIN(u,nx,ny,nz,imin,jmin,kmin)
c***
c     get the min value, and its indices, of a scalar field
c***
      INCLUDE 'undefined.h'
      INTEGER imin, jmin, kmin, nx, ny, nz, i, j, k 
      DOUBLE PRECISION u(nx,ny,nz)
c***
      FINDMIN = u(2,2,2)
      imin = 2
      jmin = 2
      kmin = 2
      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               if (u(i,j,k) .lt. FINDMIN) then
                  FINDMIN = u(i,j,k)
                  imin = i
                  jmin = j
                  kmin = k
               endif
            enddo
         enddo
      enddo
c***
      return
      end
c***********************************************************************
      SUBROUTINE finddiff(u,v,w,y,nx,ny,nz,imax,jmax,kmax)
c***
c     get the max difference, and its indices, between the components 
c     of a vector field
c***
      INCLUDE 'undefined.h'
      INTEGER imax, jmax, kmax, nx, ny, nz, i, j, k 
      DOUBLE PRECISION u(nx,ny,nz), v(nx,ny,nz), w(nx,ny,nz),
     %                 y, res, t1, t2, t3
      INTRINSIC DABS,DMAX1 
c***
      y=0.d0
      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               t1=DABS(u(i,j,k)-v(i,j,k))
               t2=DABS(u(i,j,k)-w(i,j,k))
               t3=DABS(v(i,j,k)-w(i,j,k))
               res = DMAX1(t1,t2,t3)
               if(res .gt. y) then
                  y = res
                  imax = i
                  jmax = j
                  kmax = k
               endif
            enddo
         enddo
      enddo
c***
      return
      end
c***********************************************************************
      DOUBLE PRECISION FUNCTION SUMFIELD(u,nx,ny,nz)
c***
c     get the sum of a scalar field
c***
      INCLUDE 'undefined.h'
      INTEGER nx, ny, nz, i, j, k
      DOUBLE PRECISION u(nx,ny,nz)
c***
      SUMFIELD=0.d0
      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               SUMFIELD = SUMFIELD + u(i,j,k) 
            enddo
         enddo
      enddo
c***
      return
      end
c***********************************************************************
      DOUBLE PRECISION FUNCTION FINDN2(u,v,w,nx,ny,nz)
c***
c     get the average module of a vector field
c***
      INCLUDE 'undefined.h'
      INTEGER nx, ny, nz, i, j, k
      DOUBLE PRECISION u(nx,ny,nz), v(nx,ny,nz), w(nx,ny,nz), tnpt
      INTRINSIC DSQRT
c***
      FINDN2=0.d0
      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               FINDN2 = FINDN2 + u(i,j,k)**2 + v(i,j,k)**2 + w(i,j,k)**2
            enddo
         enddo
      enddo 
      tnpt = DFLOAT((nx-2)*(ny-2)*(nz-2))
      FINDN2 = DSQRT(FINDN2/tnpt)
c***
      return
      end
c***********************************************************************
      DOUBLE PRECISION FUNCTION FINDN2S(u,nx,ny,nz)
c***
c     get the L_2 Norm of a scalar field
c***
      INCLUDE 'undefined.h'
      INTEGER nx, ny, nz, i, j, k
      DOUBLE PRECISION u(nx,ny,nz), tnpt
      INTRINSIC DSQRT
c***
      FINDN2S=0.d0
      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               FINDN2S = FINDN2S + u(i,j,k)**2 
            enddo
         enddo
      enddo 
      tnpt = DFLOAT((nx-2)*(ny-2)*(nz-2))
      FINDN2S = DSQRT(FINDN2S/tnpt)
c***
      return
      end
c***********************************************************************
      DOUBLE PRECISION FUNCTION FINDNINF(u,v,w,nx,ny,nz)
c***
c     get the max module of a vector field
c***
      INCLUDE 'undefined.h'
      INTEGER nx, ny, nz, i, j, k
      DOUBLE PRECISION u(nx,ny,nz), v(nx,ny,nz), w(nx,ny,nz), vv
      INTRINSIC DSQRT,DMAX1
c***
      FINDNINF=0.d0
      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               vv = DSQRT(u(i,j,k)**2 + v(i,j,k)**2 + w(i,j,k)**2)
               FINDNINF = DMAX1(vv,FINDNINF)
            enddo
         enddo
      enddo
c***
      return
      end
c***********************************************************************
      SUBROUTINE fill0(u,nx,ny,nz)
c***
c     fill of zeroes a 3D array
c***
      INCLUDE 'undefined.h'
      INTEGER nx,ny,nz,i,j,k
      DOUBLE PRECISION u(nx,ny,nz)
c***
      do k=1,nz
         do j=1,ny
            do i=1,nx
               u(i,j,k)=0.d0
            enddo
         enddo
      enddo
c***
      return
      end
c***********************************************************************
      SUBROUTINE copy(aout,ain,nx,ny,nz)
c***
c     copy ain in aout
c***
      INCLUDE 'undefined.h'
      INTEGER nx,ny,nz,i,j,k
      DOUBLE PRECISION ain(nx,ny,nz),aout(nx,ny,nz)
c***
      do k=1,nz
         do j=1,ny
            do i=1,nx
               aout(i,j,k)=ain(i,j,k)
            enddo
         enddo
      enddo
c***
      return
      end
c***********************************************************************
      DOUBLE PRECISION FUNCTION CFL(u,v,w,nx,ny,nz)
c***
c    check the CFL condition for the interface motion
c***
      INCLUDE 'undefined.h'
      INTEGER nx, ny, nz, i, j, k
      DOUBLE PRECISION u(nx,ny,nz), v(nx,ny,nz), w(nx,ny,nz)
      INTRINSIC DMAX1, DABS
c***
      CFL = 0.0d0
      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               CFL = DMAX1(CFL,DABS(u(i,j,k)))
            enddo
         enddo
      enddo
      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               CFL = DMAX1(CFL,DABS(v(i,j,k)))
            enddo
         enddo
      enddo
      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               CFL = DMAX1(CFL,DABS(w(i,j,k)))
            enddo
         enddo
      enddo
c***
      return
      end