program OW_editvel_mut

use sep

!implicit none

logical			:: vmut,vper,vvel,vchg
character*80                  :: hor1
real			:: o1,o2,d1,d2,v,vperc
integer			:: i1,i2,n1,n2
integer,allocatable           :: a(:)
real, allocatable	          :: pick1(:)
real,allocatable		:: velin(:,:),vsub(:,:)

call sep_init()
call sep_begin_prog()

call from_history("n1",n2); call from_history("d1",d2); call from_history("o1",o2)
call from_history("n3",n1); call from_history("d3",d1); call from_history("o3",o1)
call to_history("n1",n2); call to_history("d1",d2); call to_history("o1",o2)
call to_history("n2",1); call to_history("d2",1.); call to_history("o2",0.)
call to_history("n3",n1); call to_history("d3",d1); call to_history("o3",o1)
call to_history("n4",1); call to_history("d4",1.); call to_history("o4",0.)
call to_history("eout",4);call to_history("label3","")

allocate(velin(n2,n1));velin=0;
call sreed("in",velin,4*n1*n2)

call from_param("vmut",vmut,.false.)
call from_param("vper",vper,.false.)
call from_param("vvel",vvel,.false.) 
call from_param("vchg",vchg,.false.) 

if (vmut) call from_param("v",v)
if (vchg) call from_param("v",v)
if (vper) call from_param("vperc",vperc)
if (vvel) then
   allocate(vsub(n2,n1)) 
   call sreed("vsub",vsub,4*n1*n2)
end if

!!!!!!!
! horizon
if (.not.vchg) then
   call from_param('hor1',hor1)
   open(51,file=hor1,status='old')
   allocate(pick1(n2));pick1=9999.
   do i2=1,n2
      read(51,*,end=1) pick1(i2)
   end do
1  close(51)
end if
if (.not.vchg) then
   outer: do i2=1,n2
      do i1=1,n1
         if (i1>=pick1(i2)) then
            if (vmut) velin(i2,i1)=v
            if (vper) velin(i2,i1)=velin(i2,i1)*vperc
            if (vvel) velin(i2,i1)=vsub(i2,i1)
!            cycle outer
         end if
      end do
   end do outer
else
   allocate(a(n2));a=99999
   do i2=1,n2
      do i1=2,n1-3
         if (maxval(velin(i2,i1:i1+1))==3900.and.minval(velin(i2,i1:i1+1))==3900) then
            a(i2)=i1
            goto 10
         end if
!         if (velin(i2,i1)<velin(i2,i1-1).and.velin(i2,i1-1) > v) then
!            a(i2)=i
!            goto 10
!         else if (velin(i2,i1)==4500.and.velin(i2,i1-1)==3900) then
!            if (i1-1 < a(i2)) a(i2)=i1-1
!            goto 10
!         end if
      end do
10 end do
!   do i2=1,n2
!      do i1=2,n1
!         if (velin(i2,i1)==3900.and.a(i2)==99999) a(i2)=i1
!      end do
!   end do
   do i2=1,n2
      if (a(i2).ne.99999) then
         velin(i2,a(i2):n1)=velin(i2,a(i2))*.95
      end if
   end do
write(0,*) "passei01"
   deallocate(a)
end if
call srite("out",velin,4*n1*n2)
deallocate(velin)
if (vvel) deallocate(vsub)
call sep_end_prog()
end program
