program BinChevron
  ! Matchwell n1=100 n2=30 opertype=1     \
  ! 	binwell=binwell.H binseis=binseis.H \
  !   smseis=smseis.H misfit=misfit.H       > result.H

  use sep
  use normbin2_mod

  implicit none

  integer                               :: n12, n1, n2, nd_well, nd_seis
  integer                               :: stat
  real,    dimension (:)  , pointer     :: gridseis, gridwell, seisfold, wellfold
  real,    dimension (:)  , pointer     :: z_well, z_seis
  real,    dimension (:,:), pointer     :: xy_well, xy_seis
  real	                                :: xmin1, xmin2, d1, d2

  call sep_init(SOURCE)
  !--------------------------------------------------------------------------!
  ! Parameter Reading/Writing                                                !
  !--------------------------------------------------------------------------!
  call from_param( "nd_well", nd_well)                                     !
  call from_param( "nd_seis", nd_seis)                                     !
  call from_param( "n1", n1, 100)                                          !
  call from_param( "n2", n2, 30)                                           !
  call from_either( "xmin", xmin1)                                         !
  call from_either( "ymin", xmin2)                                         !
  call from_either( "d1", d1)                                              !
  call from_either( "d2", d2)                                              !
  !
  call to_history( "o1", xmin1)                                            !
  call to_history( "d1", d1)                                               !
  call to_history( "o2", xmin2)                                            !
  call to_history( "d2", d2)                                               !
  !
  n12 = n1 * n2                                                            !
  !--------------------------------------------------------------------------!

  !--------------------------------------------------------------------------!
  ! Array Allocation                                                         !
  !--------------------------------------------------------------------------!
  allocate( z_well(nd_well),    z_seis(nd_seis)      )                     !
  allocate( xy_well(nd_well,2), xy_seis(nd_seis,2)   )                     !
  allocate( gridseis(n12), gridwell(n12) )                                 !
  allocate( seisfold(n12), wellfold(n12) )                                 !
  !---------------------------------------------------------------------------

  !--------------------------------------------------------------------------!
  ! Read seismic and well data directly from ASCII textfile.                 !
  !--------------------------------------------------------------------------!
  open( unit=10, file='s2000.5', status='old')                             !
  do stat = 1,nd_seis                                                      !
     read( 10, 1000)  xy_seis(stat,1), xy_seis(stat,2), z_seis(stat)        !
  end do                                                                   !
1000 format( 15x, f13.6, 5x, f13.6, 6x, f12.6)                             !
  close( unit=10)                                                          !
  !
  open( unit=20, file='w2000.5', status='old')                             !
  do stat = 1,nd_well                                                      !
     read( 20, 2000)  xy_well(stat,1), xy_well(stat,2), z_well(stat)        !
  end do                                                                   !
2000 format( 15x, f13.6, 5x, f13.6, 6x, f12.6)                             !
  close( unit=20)                                                          !
  !
  call sep_close ()                                                        !
  !--------------------------------------------------------------------------!

  !--------------------------------------------------------------------------!
  ! By convention, set depths as positive numbers                            !
  !--------------------------------------------------------------------------!
  z_well = -z_well;  z_seis = -z_seis                                      !
  !--------------------------------------------------------------------------!

  !--------------------------------------------------------------------------!
  ! Convert (x,y) to km from feet.                                           !
  !--------------------------------------------------------------------------!
  xy_well = xy_well/3280.1;  xy_seis = xy_seis/3280.1;                     !
  !--------------------------------------------------------------------------!

  !--------------------------------------------------------------------------!
  ! Interpolate (nearest-neighbor binning) well and seismic data to a        !
  ! regular grid, also normalizing by the fold.                              !
  !--------------------------------------------------------------------------!
  call normbin2_init( n1, n2, nd_well, xmin1, xmin2, d1, d2, xy_well )     !
  stat = normbin2_lop( .true., .false., gridwell, z_well )                 !
  call get_fold( wellfold )                                                !
  !
  call normbin2_init( n1, n2, nd_seis, xmin1, xmin2, d1, d2, xy_seis )     !
  stat = normbin2_lop( .true., .false., gridseis, z_seis )                 !
  call get_fold( seisfold )                                                !
  !--------------------------------------------------------------------------!

  !--------------------------------------------------------------------------!
  ! Write binned seismic and well data (fold on frame 2) to disk             !
  !--------------------------------------------------------------------------!
  call to_history( "n1", n1,    "binwell" )                                !
  call to_history( "n2", n2,    "binwell" )                                !
  call to_history( "n3", 2,     "binwell" )                                !
  call to_history( "o1", xmin1, "binwell" )                                !
  call to_history( "o2", xmin2, "binwell" )                                !
  call to_history( "d1", d1,    "binwell" )                                !
  call to_history( "d2", d2,    "binwell" )                                !
  call sep_write( gridwell, "binwell")                                     !
  call sep_write( wellfold, "binwell")                                     !
  call to_history( "meanval", sum(z_well)/nd_well, "binwell" )             !
  call to_history( "nd_well", nd_well, "binwell" )                         !
  !
  call to_history( "n1", n1,    "binseis" )                                !
  call to_history( "n2", n2,    "binseis" )                                !
  call to_history( "n3", 2,     "binseis" )                                !
  call to_history( "o1", xmin1, "binseis" )                                !
  call to_history( "o2", xmin2, "binseis" )                                !
  call to_history( "d1", d1,    "binseis" )                                !
  call to_history( "d2", d2,    "binseis" )                                !
  call sep_write( gridseis, "binseis")                                     !
  call sep_write( seisfold, "binseis")                                     !
  call to_history( "meanval", sum(z_seis)/nd_seis, "binseis" )             !
  call to_history( "nd_seis", nd_seis, "binseis" )                         !
  !--------------------------------------------------------------------------!
  if( exist_file("xyzseis") ) then
     call to_history( "n1", nd_seis, "xyzseis" )                            !
     call to_history( "n2", 3,       "xyzseis" )                            !
     call sep_write( xy_seis, "xyzseis" )
     call sep_write( z_seis, "xyzseis" )
  end if
  if( exist_file("xyzwell") ) then
     call to_history( "n1", nd_well, "xyzwell" )                            !
     call to_history( "n2", 3,       "xyzwell" )                            !
     call sep_write( xy_well, "xyzwell" )
     call sep_write( z_well, "xyzwell" )
  end if
  !--------------------------------------------------------------------------!

  call exit(0)

end program BinChevron

