program Model2d
  !! Create a 2.5-D velocity model
  !! Paul Sava (paul@sep.stanford.edu)
  !! 09/19/00

  use sep

  implicit none
  character (len=10) :: buffer

  integer :: nz,  nx,  ny
  real    :: oz,  ox,  oy
  real    :: dz,  dx,  dy
  real    :: zmin,xmin
  real    :: zmax,xmax
  real    :: z,   x

  real, dimension(:,:), allocatable:: model

  integer :: ne !! elements in the model
  integer :: ie

  real    :: vback

  real, dimension(:), allocatable :: xs,zs !! start coords
  real, dimension(:), allocatable :: xe,ze !! end   coords
  real, dimension(:), allocatable :: v !! velocity
  real, dimension(:), allocatable :: r !! reflectivity

  integer :: ixs,izs, ix,jx
  integer :: ixe,ize, iz,jz
  integer :: iy

  call sep_init()

  call to_history("esize",4)

  call from_param("nz",nz);   call to_history("n1",nz)
  call from_param("oz",oz,0.);call to_history("o1",oz)
  call from_param("dz",dz,.1);call to_history("d1",dz)
  zmin = oz
  zmax = oz+(nz-1)*dz

  call from_param("nx",nx);   call to_history("n2",nx)
  call from_param("ox",ox,0.);call to_history("o2",ox)
  call from_param("dx",dx,.1);call to_history("d2",dx)
  xmin = ox
  xmax = ox+(nx-1)*dx

  call from_param("ny",ny,1); call to_history("n3",ny)
  call from_param("oy",oy,0.);call to_history("o3",oy)
  call from_param("dy",dy,1.);call to_history("d3",dy)

  call from_param("vback",vback,0.0)
  call from_param("ne",ne,1)

  allocate(model(nz,nx))
  model=vback !! init everything with the background

  allocate(xs(ne), zs(ne))
  allocate(xe(ne), ze(ne))
  allocate(v(ne))
  allocate(r(ne))

  do ie=1,ne
     !! X axis
     if(ie<10) then
        write (buffer, "(a,i1)" ) 'xs',ie !! start
     else if(ie<100) then
        write (buffer, "(a,i2)" ) 'xs',ie
     else if(ie<1000) then
        write (buffer, "(a,i3)" ) 'xs',ie
     end if
     call from_param(buffer,xs(ie),xmin)
     if(xs(ie)<xmin) xs(ie)=xmin
     if(xs(ie)>xmax) xs(ie)=xmax

     if(ie<10) then
        write (buffer, "(a,i1)" ) 'xe',ie
     else if(ie<100) then
        write (buffer, "(a,i2)" ) 'xe',ie
     else if(ie<1000) then
        write (buffer, "(a,i3)" ) 'xe',ie
     end if
     call from_param(buffer,xe(ie),xmax)
     if(xe(ie)<xmin) xe(ie)=xmin
     if(xe(ie)>xmax) xe(ie)=xmax

     !! Z axis
     if(ie<10) then
        write (buffer, "(a,i1)" ) 'zs',ie !! start
     else if(ie<100) then
        write (buffer, "(a,i2)" ) 'zs',ie
     else if(ie<1000) then
        write (buffer, "(a,i3)" ) 'zs',ie
     end if
     call from_param(buffer,zs(ie),zmin)
!     if(zs(ie)<zmin) zs(ie)=zmin
!     if(zs(ie)>zmax) zs(ie)=zmax

     if(ie<10) then
        write (buffer, "(a,i1)" ) 'ze',ie
     else if(ie<100) then
        write (buffer, "(a,i2)" ) 'ze',ie
     else if(ie<1000) then
        write (buffer, "(a,i3)" ) 'ze',ie
     end if
     call from_param(buffer,ze(ie),zmax)
!     if(ze(ie)<zmin) ze(ie)=zmin
!     if(ze(ie)>zmax) ze(ie)=zmax

     if(zs(ie)<zmin .or. zs(ie)>zmax) cycle
     if(ze(ie)<zmin .or. ze(ie)>zmax) cycle

     !! velocity
     if(ie<10) then
        write (buffer, "(a,i1)" ) 'v',ie
     else if(ie<100) then
        write (buffer, "(a,i2)" ) 'v',ie
     else if(ie<1000) then
        write (buffer, "(a,i3)" ) 'v',ie
     end if
     call from_param(buffer,v(ie),vback)
     !! reflectivity
     if(ie<10) then
        write (buffer, "(a,i1)" ) 'r',ie
     else if(ie<100) then
        write (buffer, "(a,i2)" ) 'r',ie
     else if(ie<1000) then
        write (buffer, "(a,i3)" ) 'r',ie
     end if
     call from_param(buffer,r(ie),0.0)

     ixs = 1+nint((xs(ie)-xmin)/dx)
     ixe = 1+nint((xe(ie)-xmin)/dx)

     izs = 1+nint((zs(ie)-zmin)/dz)
     ize = 1+nint((ze(ie)-zmin)/dz)

     if(xe(ie)-xs(ie) /=0.) then
        if(r(ie)/=0) then !! reflectivity
           do jx=ixs,ixe              
              x = xmin + (jx-1)*dx

              z = zs(ie)
              z = z + (ze(ie)-zs(ie))/(xe(ie)-xs(ie)) * (x-xs(ie))
              jz = 1+ nint((z-zmin)/dz)
              model(jz,jx) = r(ie)
           end do
        else
           do jx=ixs,ixe
              x = xmin + (jx-1)*dx

              z = zs(ie)
              z = z + (ze(ie)-zs(ie))/(xe(ie)-xs(ie)) * (x-xs(ie))

              iz = 1+ nint((z-zmin)/dz)
              do jz=iz,nz
                 model(jz,jx) = v(ie)
              end do
           end do
        end if
     else
        if(r(ie)/=0) then !! reflectivity
           do jz=izs,ize
              model(jz,ixs) = r(ie)
           end do
        else
           do jz=izs,nz
              model(jz,ixs) = v(ie)
           end do
        end if
     end if



  end do

  do iy=1,ny
     call sep_write(model)
  end do
  deallocate(model)

  call sep_close()
  call exit(0)
end program MODEL2D













