c   $Id: printvel.f,v 1.1 1998/10/20 13:23:03 zaleski Exp $
c   (c) S. Zaleski, CNRS, UPMC
c
      SUBROUTINE initpvel(nx,ny)
      include 'undefined.h'
      integer nx,ny,i,system,record
      include 'machdep.f'
      common /pvel/record
      include 'mach.h'
      record=1
C
      if(cray.eq.0) then
      OPEN(unit=14,form="unformatted",file="velocity",access="direct",
     *recl=nx*ny*16/reclunit)
      close(unit=14,status="delete")
      OPEN(unit=14,form="unformatted",file="velocity",access="direct",
     *status="new",recl=nx*ny*16/reclunit)
      OPEN(unit=16,file="front")
      close(unit=16,status="delete") 
      OPEN(unit=16,file="front",status="new")
      OPEN(unit=18,file="profile")
      close(unit=18,status="delete") 
      OPEN(unit=18,file="profile",status="new")
      endif
      RETURN      
      END

      SUBROUTINE printvel(tau,u,v,c,p,su,sv,sc,sp,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny,record,i,j
      DOUBLE PRECISION u(nx,ny),v(nx,ny),c(nx,ny),p(nx,ny)
      DOUBLE PRECISION tau,tauh
      REAL su(nx,ny),sv(nx,ny),sc(nx,ny),sp(nx,ny)
      common /pvel/record

      include 'mach.h'
      if(cray.eq.0) then
      do j=1,ny
         do i=1,nx
            su(i,j) = SNGL(u(i,j))
            sv(i,j) = SNGL(v(i,j))
            sc(i,j) = SNGL(c(i,j))
 	    sp(i,j) = SNGL(p(i,j))
        enddo
      enddo
      tauh=1.0d0/(tau*(nx-2))
      do j=1,ny
         do i=1,nx
            su(i,j)=SNGL(su(i,j)*tauh)
            sv(i,j)=SNGL(sv(i,j)*tauh)
            sp(i,j)=SNGL(sp(i,j)*tauh**2)
         enddo
      enddo
c
      write(14,rec=record) su,sv,sc,sp
      record = record + 1
      call flush(14)
      endif
      do i=1,nx
         write(18,*) i,u(i,2)
      enddo
      write(18,*) ' '
      return
      end


      SUBROUTINE mwrit(c,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny,i,j
      DOUBLE PRECISION c(nx,ny)
      do j=1,ny
         do i=1,nx
            write (54,*) c(i,j)
         enddo
      enddo
      return
      end

      SUBROUTINE front(c,nx,ny,tphys)
      include 'undefined.h'
      INTEGER nx,ny,compt,i,j
      DOUBLE PRECISION c(nx,ny),h,eps,tphys
      
      eps=1E-6
      h=1.d0/(nx-2)
         do i=1,nx-1
            if(c(nx-i,2).gt.0.8) then
               write(16,*) tphys,h*(nx-i-2)
               go to 2
         endif
      enddo
 2    continue
      return
      end