program OW_editvel_mud

use sep

implicit none
real			:: o1,o2,d1,d2,vh2o,vvv
integer			:: i1,i2,n1,n2
integer, allocatable	:: 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),ind1(n2));ind1=0;velin=0;
call sreed("in",velin,4*n1*n2)

call from_param("vh2o",vh2o)
open(unit=2,file='dump')
do i2=1,n2
   do i1=1,300
      if (int(velin(i1,i2)).eq.int(vh2o)) then
         ind1(i2)=i1
      else
         exit
      end if
   end do
end do
if (exist_file('wbott')) then
   call to_history("n1",n2,'wbott'); call to_history("d1",d2,'wbott'); call to_history("o1",o2,'wbott')
   call to_history("n2",1,'wbott'); call to_history("d2",1.,'wbott'); call to_history("o2",1.,'wbott')
   call to_history("eout",4,'wbott');
   call srite('wbott',float(ind1)*d1,4*n2)
end if
do i2=1,n2
   do i1=1,n1
      if (i1 <= ind1(i2)) then
         velin(i1,i2) = vh2o
      else if (i1 > ind1(i2) .and. i1 <= ind1(i2)+int(405/d1) ) then
         velin(i1,i2) = 1580
      else if (i1 > ind1(i2)+int(405/d1) .and. i1 <= ind1(i2)+int(1120/d1)) then
         velin(i1,i2) = 0.5*(i1-ind1(i2)-int(405/d1)-1)*d1+1580
      else if (i1 > ind1(i2)+int(1120/d1) .and. i1 <= ind1(i2)+int(2135/d1)) then
         velin(i1,i2) = 1.0*(i1-ind1(i2)-int(1120/d1)-1)*d1+velin(ind1(i2)+int(1120/d1),i2)
      else if (i1 > ind1(i2)+int(2135/d1) .and. i1 <= ind1(i2)+int(4035/d1)) then
         velin(i1,i2) = 0.6*(i1-ind1(i2)-int(2135/d1)-1)*d1+velin(ind1(i2)+int(2135/d1),i2)
      else if (i1 > ind1(i2)+int(4035/d1) .and. i1 <= n1) then
         velin(i1,i2) = 0.1*(i1-ind1(i2)-int(4035/d1)-1)*d1+velin(ind1(i2)+int(4035/d1),i2)
      end if
   end do
end do
call srite("out",velin,4*n1*n2)
deallocate(velin)
close(2)
call sep_end_prog()
end program
