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
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***********************************************************************
      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***********************************************************************
      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 BOXSUMFIELD(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***
      BOXSUMFIELD=0.d0
      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               BOXSUMFIELD = BOXSUMFIELD + u(i,j,k) 
            enddo
         enddo
      enddo
c***
      return
      end