program OW_editvel_mud1

use sep

implicit none
real			:: o1,o2,d1,d2,om1,om2,dm1,dm2,vvv,r,grad
integer			:: i1,i2,n1,n2,nm1,nm2,wn1,wn2,nzero,ia,ib
real,allocatable		:: velin(:,:),mask1(:,:),vsub(:),vtemp(:),mask2(:,:)

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","")

call from_aux("mask1","n1",nm1); call from_aux("mask1","d1",dm1); call from_aux("mask1","o1",om1)
call from_aux("mask1","n2",nm2); call from_aux("mask1","d2",dm2); call from_aux("mask1","o2",om2)

allocate(velin(n1,n2));velin=0;
allocate(mask1(nm1,nm2));mask1=0;
allocate(mask2(nm1,nm2));mask2=0;
allocate(vtemp(nm2));vtemp=0;

call sreed("in",velin,4*n1*n2)
call sreed("mask1",mask1,4*nm1*nm2)
open(unit=2,file='dump')

do i1=1,nm1
   where(mask1(i1,:).ge.50.)
      mask1(i1,:)=1
   elsewhere
      mask1(i1,:)=0
   end where
end do
!r=1.05
!mask2=mask1*r
!do i1=1,nm1
!   where(mask2(i1,:).eq.0.)
!      mask2(i1,:)=1
!   end where
!end do
!
!do i2=1,nm2
!do i1=1,nm1
!   if (mask2(i1,i2).eq.r.and.i2.lt.160) then
!      mask2(i1,i2)=0.92
!   else if (mask2(i1,i2).eq.r.and.i2.gt.160.and.i2.lt.240) then
!      mask2(i1,i2)=0.95
!   else if (mask2(i1,i2).eq.r.and.i2.gt.240.and.i2.lt.360) then
!      mask2(i1,i2)=1.01
!   else if (mask2(i1,i2).eq.r.and.i2.gt.490) then
!      mask2(i1,i2)=1.11
!   end if
!end do
!end do

call to_history("n1",nm1,"maskout"); call to_history("d1",dm1,"maskout"); call to_history("o1",om1,"maskout")
call to_history("n2",nm2,"maskout"); call to_history("d2",dm2,"maskout"); call to_history("o2",om2,"maskout")
call srite("maskout",mask1,4*nm1*nm2)

wn1=int((om1-o1)/d1)+1
wn2=int((om2-o2)/d2)+1

!do i1=1,nm1
!   vtemp=velin(i1+wn1,wn2:wn2+nm2-1)
!   nzero = count(mask = (mask1(i1,:).ne.1))
!   allocate(vsub(nzero))
!   where(mask1(i1,:).ne.1)
!      vsub=vtemp
!   end where
!   vvv=sum(vsub)/nzero
!   where(mask1(i1,:).eq.1)
!      vtemp=vvv
!   end where
!   vtemp=vtemp*mask2(i1,:)
!   velin(i1+wn1,wn2:wn2+nm2-1)=vtemp
!   deallocate(vsub)
!end do


do i2=1,nm2
   ia=0;ib=1
   do i1=1,nm1
      if (mask1(i1,i2)==1.and.ia==0) then
         ia=i1-15
         ib=ia
      end if
      if (mask1(i1,i2)==0.and.ib==ia) then
         ib=i1+15
      end if
   end do
   write(2,*) i2,ia,ib
   if (sum(mask1(:,i2)) > 0) then
      grad=(velin(ib+wn1,i2+wn2)-velin(ia+wn1,i2+wn2))/(ib-ia)
      write(0,*) sum(mask1(:,i2)),i2,grad
   end if
   do i1=1,nm1
      if (i1.ge.ia.and.i1.le.ib) then
         velin(i1+wn1,i2+wn2)=grad*(i1-ia)+velin(ia+wn1,i2+wn2)
      end if
   end do
end do

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