 program mute1

 use sep

 implicit none

  type axis
     integer 					:: n
     real    					:: o, d, v ! v=(i-1)*d+o
     real, allocatable				:: all(:)
  end type axis

  real, allocatable				:: data(:,:,:),vel(:,:),mask(:,:),maxvel(:)
  integer					:: ih,ix,iz
  real						:: cx,cz,maxw
  type(axis)					:: x,z,h,vx,vz

 call sep_init()
 call sep_begin_prog()

   call from_aux("cig","n1",z%n); call from_aux("cig","d1",z%d); call from_aux("cig","o1",z%o)
   call from_aux("cig","n2",h%n); call from_aux("cig","d2",h%d); call from_aux("cig","o2",h%o)
   call from_aux("cig","n3",x%n); call from_aux("cig","d3",x%d); call from_aux("cig","o3",x%o)

   call from_aux("vel","n1",vz%n); call from_aux("vel","d1",vz%d); call from_aux("vel","o1",vz%o)
   call from_aux("vel","n2",vx%n); call from_aux("vel","d2",vx%d); call from_aux("vel","o2",vx%o)

   call to_history("n1",z%n); call to_history("d1",z%d); call to_history("o1",z%o)
   call to_history("n2",h%n); call to_history("d2",h%d); call to_history("o2",h%o)
   call to_history("n3",x%n); call to_history("d3",x%d); call to_history("o3",x%o)

   if (vz%n.ne.z%n.or.vx%n.ne.x%n) call seperr('Different dimensions: fix it')
   allocate(vel(vz%n,vx%n),mask(vz%n,vx%n),data(z%n,h%n,x%n),maxvel(z%n))
   call sreed("vel",vel, 4*vz%n*vx%n)

   mask=0;where(vel.lt.14000) mask=1

   do ix=1,x%n
      cx=(ix-1)*x%d+x%o
      if (cx < 34700) then
	 mask(:,ix)=0.
      else
         do iz=z%n,1,-1
            if (mask(iz,ix) .eq. 0.) then
	       mask(iz:iz+10,ix-10:ix)=0.
    	       cycle
	    end if
         end do
      end if
   end do

   do iz=1,z%n
      cz=(iz-1)*z%d+z%o
      if (cz < 9500) mask(iz,:)=0.
   end do

   do ix=1,x%n; do ih=1,h%n; do iz=1,z%n
      data(iz,ih,ix)=mask(iz,ix)
   end do;end do;end do


   call srite("out",data, 4*z%n*h%n*x%n)
   deallocate(data)

 call sep_end_prog()
 end program mute1
