program OW_vrms
!
! dim1 = depth
! dim2 = xline
! dim3 = iline
!
use sep

implicit none
logical                       :: int2rms
real			:: o1,o2,o3,d1,d2,d3
integer			:: i1,i2,i3,n1,n2,n3
real,allocatable		:: velin(:,:),velout(:,:)

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 from_param("int2rms",int2rms,.true.)
allocate(velin(n1,n2));velin=0;
allocate(velout(n1,n2));velout=0;

if(int2rms) then
  do i3=1,n3
     call sreed("in",velin,4*n1*n2)
     do i2=1,n2
     do i1=1,n1
        velout(i1,i2)=sum(velin(1:i1,i2))/i1
     end do
     end do
     call srite("out",velout,4*n1*n2)
  end do
else
  do i3=1,n3
     call sreed("in",velin,4*n1*n2)
     velout(1,:)=velin(1,:)
     do i2=1,n2
     do i1=2,n1
!        velout(i1,i2)=sqrt((i1*velin(i1,i2)**2-sum(velout(1:i1-1,i2)**2)))
        velout(i1,i2)=i1*velin(i1,i2)-sum(velout(1:i1-1,i2))
     end do
     end do
     call srite("out",velout,4*n1*n2)
  end do
end if
deallocate(velin,velout)
call sep_end_prog()
end program
