program OW_3deditvel_mut

use sep

!implicit none

logical			:: vmut,vper,vvel,vthr,vextr
character*80                  :: hor1
real			:: o1,o2,o3,d1,d2,d3,v,vperc,s1,dz
integer			:: i1,i2,i3,n1,n2,n3,m1
real, allocatable	          :: pick1(:)
real,allocatable		:: velin(:,:),vsub(:,:)

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 from_history("n3",n3); call from_history("d3",d3); call from_history("o3",o3)
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",n3); call to_history("d3",d3); call to_history("o3",o3)
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(n1,n2));velin=0;

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

if (vmut) call from_param("vextr",vextr,.false.)
if (not(vextr.or.vvel).or.vthr) call from_param("v",v)
if (vextr) call from_param("dz",dz,-4*d1)
if (vper) call from_param("vperc",vperc)
if (vvel) then
   allocate(vsub(n1,n2));vsub=0.
end if

!!!!!!!
! horizon
   call from_param('hor1',hor1)
   open(51,file=hor1,status='old')
   allocate(pick1(n2*n3));pick1=99999.
   do i2=1,n2*n3
      read(51,*,end=1) pick1(i2)
   end do
1  close(51)

   do i3=1,n3
      call sreed("in",velin,4*n1*n2)
      if (vvel) then
         call sreed("vsub",vsub,4*n1*n2)
      end if
      do i2=1,n2
         do i1=1,n1
            s1=(i1-1)*d1+o1
            if (s1 >= pick1((i3-1)*n2+i2)+dz) then
               if (vmut) then
                  if (vextr) then 
                     m1=((pick1((i3-1)*n2+i2)+dz)/d1)+1.5
                     velin(i1,i2)=velin(m1,i2)
                  else
                     velin(i1,i2)=v
                  end if
               end if
               if (vper) velin(i1,i2)=velin(i1,i2)*vperc
               if (vvel) velin(i1,i2)=vsub(i1,i2)
               if (vthr) then
                  if (velin(i1,i2) < v) velin(i1,i2) = v
               end if
            end if
         end do
      end do
      call srite("out",velin,4*n1*n2)
   end do

   deallocate(velin)
   if (vvel) deallocate(vsub)
   call sep_end_prog()
end program
