 program OW_mutvel

 use sep

 implicit none

 real				:: oz,ox,dz,dx,x,z,f,g
 real,allocatable			:: vel(:,:),vec1(:),vec2(:),vec3(:)
 integer				:: ix, iz, ih, nz,nx,fi,dir

 call sep_init()
 call sep_begin_prog()

 call from_history("n1",nz);call from_history("d1",dz);call from_history("o1",oz)
 call from_history("n2",nx);call from_history("d2",dx);call from_history("o2",ox)

 call to_history("n1",nz);call to_history("d1",dz);call to_history("o1",oz)
 call to_history("n2",nx);call to_history("d2",dx);call to_history("o2",ox)

 call from_param("dir",dir)

 allocate(vel(nz,nx))

 call sreed("in",vel,4*nz*nx)

 if (dir==2) then
    allocate(vec1(nz),vec2(nz),vec3(nz))
    vec1=vel(:,int((14000-ox)/dx)+1)
    vec2=vel(:,int((44000-ox)/dx)+1)
    do ix=1,nx
       x = (ix-1)*dx+ox
       if (x>14000 .and. x<44000) then
          f=(x-14000)/30000;fi=f;f=f-fi;g=1-f
          vec3=g*vec1+f*vec2
          vel(:,ix)=vec3
       end if
    end do
else
    allocate(vec1(nx),vec2(nx),vec3(nx))
    vec1=vel(int((900-oz)/dz)+1,:)
    vec2=vel(int((1600-oz)/dz)+1,:)
    do iz=1,nz
       z = (iz-1)*dz+oz
       if (z>900 .and. z<1600) then
          f=(z-900)/1600;fi=f;f=f-fi;g=1-f
          vec3=g*vec1+f*vec2
          vel(iz,:)=vec3
       end if
    end do
end if
 call srite("out",vel,4*nz*nx)
 deallocate(vel,vec1,vec2,vec3)

 call sep_end_prog()

 end program
