program OW_editvel

use sep
use lint

implicit none
real			:: gradz,v0,vsalt,o1,o2,d1,d2,vsal,vh2o
integer			:: i1,i2,n1,n2,j
integer, allocatable	:: ind(:),ind1(:)
real,allocatable		:: velin(:,:)

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("eout",4);call to_history("label3","")

allocate(velin(n1,n2),ind(n2),ind1(n2))
call sreed("in",velin,4*n1*n2)

call from_param("gradz",gradz)
call from_param("v0",v0)
call from_param("vsalt",vsalt)
call from_param("vh2o",vh2o)

do i2=1,n2
   do i1=1,n1
      if (velin(i1,i2).eq.vsalt.and.i1.lt.1185) then
         ind(i2)=i1
      end if
      if (velin(i1,i2).le.vh2o) then
         ind1(i2)=i1
      end if
   end do
end do

do i2=1,n2
   do i1=1,n1
      if ((i1.ge.ind(i2)).and.(ind(i2).ne.0)) then
         velin(i1,i2)=gradz*(i1-1)+v0
      else if (ind(i2).eq.0.and.i1.ge.ind1(i2)+240) then
!         velin(i1,i2)=gradz*(i1-ind1(i2)-50)+1.6*v0
         velin(i1,i2)=gradz*(i1-ind1(i2)-600)+3*v0
      end if
   end do
end do

call srite("out",velin,4*n1*n2)
deallocate(velin)
call sep_end_prog()
end program

!character*80		:: infile
!real			:: vec(2)
!real,pointer		:: t(:),weight(:)
!real			:: xcoord(200),zcoord(200)
!real,allocatable,dimension(:) :: zout,vtemp
!integer			:: i,j,stat,n2
!write(0,*) "Input the horizon pick file"
!read(0,*) infile
!open(unit=2,file=infile)
!xcoord=0.;zcoord=0.; vec=0.; j=2
!do i=1,9999
!   read(2,*,end=10) vec 
!   zcoord(j)=real(vec(1)); xcoord(j)=real(vec(2)); j=j+1
!write(0,*) vec
!end do
!10 continue
!xcoord(1)=1; xcoord(j)=real(n2)
!zcoord(1)=zcoord(2); zcoord(j)=zcoord(j-1)
!allocate(weight(j),t(j),zout(n2)); weight=1.; zout=0.
!t=xcoord(1:j)
!write(0,*) "passei 00"
!call lint_init (.true.,1.,1.,n2,t,weight)
!write(0,*) "passei 01"
!stat = lint_lop (.true.,.false., zout, zcoord(1:j))
!write(0,*) "passei 02"
!zout(xcoord(j-1):n2)=zcoord(j-1)
!call srite("out",zout,4*n2)
!deallocate(weight,t,zout,vtemp)

   
