program GrMa_
  !! Compute the gradient magnitude of a velocity model
  !! Paul Sava (paul@sep.stanford.edu)
  !! 11/10/03

  !! input:  z-x 
  !! output: z-x 

  use sep
  implicit none

  real, dimension(:,:), allocatable :: data
  real, dimension(:,:), allocatable :: grad

  integer :: nz,  nx
  integer :: iz,  ix
  real    :: oz,  ox
  real    :: dz,  dx

  call sep_init()

  call from_history("n1",nz)
  call from_history("o1",oz,0.)
  call from_history("d1",dz,1.)

  call from_history("n2",nx)
  call from_history("o2",ox,0.)
  call from_history("d2",dx,1.)

  call sep_close()

  allocate(data(nz,nx)); data=0
  allocate(grad(nz,nx)); grad=0.

  call sep_read(data)

  do iz=1,nz-1
     do ix=1,nx-1
!write(0,*) iz,ix
        grad(iz,ix) = sqrt ( &
        (data(iz  ,ix+1) - data(iz,ix))**2 + &
        (data(iz+1,ix  ) - data(iz,ix))**2)
     end do
  end do

  grad(:,nx)=grad(:,nx-1)
  grad(nz,:)=grad(nz-1,:)

  call sep_write(grad)

  call exit(0)
end program GRMA_













