program main
  use sep
  implicit none

  real, dimension(:,:,:), allocatable :: vold, vnew

  integer                            :: n1, n2, n3
  integer                            :: i1, i2, i3
  real                               :: o1, o2, o3
  real                               :: d1, d2, d3
  real                               :: x1, x2, x3

  call initpar()

  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)

  allocate(vold(n1,n2,n3))
  allocate(vnew(n1,n2,n3))

  call sep_read(vold)

  vnew = vold

  do i3=1, n3
     x3 = o3+(i3-1)*d3
     if (x3>11250.) then
        do i2=1, n2; x2 = o2+(i2-1)*d2
        do i1=1, n1; x1 = o1+(i1-1)*d1
           if (vold(i1,i2,i3)<14800) then
              if (x1<56750.) then
                 vnew(i1,i2,i3) = vold(1,i2,i3)
              else if (x1>=56750. .and. x1<69350 .and. x3>=16150) then
                 vnew(i1,i2,i3) = vold(1,i2,i3)
              else if (x3>=18700.) then
                 vnew(i1,i2,i3) = vold(1,i2,i3)
              end if
           end if
        end do
        end do
     end if
  end do


  call sep_write(vnew)

end program
