! module OW_dips
 program OW_dips
 use sep 

 implicit none
! contains

! subroutine owdips(vel,dips)

integer		:: n1,n2
real		:: d1,d2,o1,o2
real,allocatable	:: vel(:,:),dips(:,:)

integer			:: i,m,n,dev
!real			:: vel(:,:),dips(:,:)
real,allocatable		:: gx(:,:),gz(:,:)

 call sep_init()
 call sep_begin_prog()

 call from_history("n1",n1);	call from_history("d1",d1);	call from_history("o1",o1)
 call from_history("n2",n2);	call from_history("d2",d2);	call from_history("o2",o2)

 call to_history("n1",n1);	call to_history("d1",d1);	call to_history("o1",o1)
 call to_history("n2",n2);	call to_history("d2",d2);	call to_history("o2",o2)
 call to_history("n3",1);	call to_history("d3",1.);	call to_history("o3",0.)  

m=n1;n=n2
!m=size(vel,1)
!n=size(vel,2)

allocate(vel(m,n),dips(m,n))
allocate(gx(m,n),gz(m,n))
gx=0;gz=0

call sreed("in",vel,4*m*n)

dev=1
if (dev==1) then
   do i=1,m
      if (i==1) then
         gz(1,:)=(vel(2,:)-vel(1,:))/d1;
      elseif (i==2) then
         gz(2,:)=(vel(3,:)-vel(1,:))/(2*d1);
      elseif (i==m-1) then
         gz(i,:)=(vel(m,:)-vel(m-2,:))/d1;
      elseif (i==m) then
         gz(m,:)=(vel(m,:)-vel(m-1,:))/d1;
      else; 
         gz(i,:)=(vel(i-2,:)-8*vel(i-1,:)+8*vel(i+1,:)-vel(i+2,:))/(12*d1);
      end if
   end do
   do i=1,n
      if (i==1) then
         gx(:,1)=(vel(:,2)-vel(:,1))/d2;
      elseif (i==2) then
         gx(:,2)=(vel(:,3)-vel(:,1))/(2*d2);
      elseif (i==n-1) then
         gx(:,i)=(vel(:,n)-vel(:,n-2))/(2*d2);
      elseif (i==n) then
         gx(:,n)=(vel(:,n)-vel(:,n-1))/d2;
      else
         gx(:,i)=(vel(:,i-2)-8*vel(:,i-1)+8*vel(:,i+1)-vel(:,i+2))/(12*d2);
      end if
   end do
else
   do i=1,m
      if (i==1) then
         gz(1,:)=(vel(2,:)-vel(1,:))/d1;
      elseif (i==m) then
         gz(m,:)=(vel(m,:)-vel(m-1,:))/d1;
      else; 
         gz(i,:)=(vel(i+1,:)-vel(i,:))/d1;
      end if
   end do
   do i=1,n
      if (i==1) then
         gx(:,1)=(vel(:,2)-vel(:,1))/d2;
      elseif (i==n) then
         gx(:,n)=(vel(:,n)-vel(:,n-1))/d2;
      else
         gx(:,i)=(vel(:,i+1)-vel(:,i))/d2;
      end if
   end do
end if

dips=atan(gx*gz/(gz*gz+0.00001))*180/3.14159

call srite("out",dips,4*m*n)
call sep_end_prog()
!end souroutine
!end module
end program
