program main
  use base_struct_mod
  use base_ioutil_mod
  use hes_fileio_mod
  use hes_gather_mod
  implicit none

  type(hes_grid)                      :: hgrid
  type(axis)                          :: aphs
  real, dimension(:,:,:), allocatable :: hess, ahes
  real, dimension(:,:,:), allocatable :: vmod
  real, dimension(:)    , allocatable :: dip_agl
  real, dimension(:)    , allocatable :: sca_agl
  real, dimension(:,:,:), allocatable :: sca_map
  real                                :: vref, ws, pct
  real                                :: sca_agl_min, sca_agl_max, dsca_agl
  real                                :: phs_min, phs_max, dphs
  integer                             :: ntpx, ntpy, ntpz, nkx, nky, nkz
  integer                             :: ix_modl, iy_modl, iz_modl, stat
  logical                             :: radian, norm
  character(len=256)                  :: hess_, vmod_, idip_, isca_, msca_, test_

  real, parameter                     :: PI = 3.14159265359

  call sep_init(SOURCE)
  call getch_add_string("head=/dev/null")
  call getch_add_string("noheader=y")
  call set_no_putch()
  call sep_begin_prog()

  hess_ = "hess"
  vmod_ = "vmod"
  idip_ = "idip"
  isca_ = "isca"
  msca_ = "msca"

  test_ = "test"

  call auxin (hess_)
  call auxout(idip_)
  call auxout(isca_)
  call auxout(msca_)
  call auxout(test_)

  call hinaxis(hess_, hgrid%ax_offd, 1)
  call hinaxis(hess_, hgrid%ay_offd, 2)
  call hinaxis(hess_, hgrid%az_offd, 3)
  call hinaxis(hess_, hgrid%ax_modl, 4)
  call hinaxis(hess_, hgrid%ay_modl, 5)
  call hinaxis(hess_, hgrid%az_modl, 6)

  call from_param("ws"    , ws    ,    10.) 
  call from_param("radian", radian, .true.)
  call from_param("pct"   , pct   ,    0.5)
  call from_param("norm"  , norm  , .true.)

  if (.not.radian) ws = ws * 2 * PI

  call from_param("sca_agl_min", sca_agl_min,  0.)
  call from_param("sca_agl_max", sca_agl_max, 60.)
  call from_param("dsca_agl"   , dsca_agl   ,  1.)

  call from_param("odip_agl", hgrid%adip_agl%o,  0.)
  call from_param("ddip_agl", hgrid%adip_agl%d,  1.)
  call from_param("ndip_agl", hgrid%adip_agl%n, 360)

  call from_param("phs_min", phs_min,   0.)
  call from_param("phs_max", phs_max, 180.)
  call from_param("dphs"   , dphs   ,  0.5)

  aphs%o = phs_min
  aphs%d = dphs
  aphs%n = 1.5 + (phs_max-phs_min)/dphs

  hgrid%asca_agl%o = sca_agl_min
  hgrid%asca_agl%d = dsca_agl
  hgrid%asca_agl%n = 1.5 + (sca_agl_max-sca_agl_min)/dsca_agl

  call houaxis(isca_, hgrid%asca_agl, 1)
  call houaxis(isca_,  hgrid%ax_modl, 2)
  call houaxis(isca_,  hgrid%ay_modl, 3)
  call houaxis(isca_,  hgrid%az_modl, 4)

  call houaxis(idip_, hgrid%adip_agl, 1)
  call houaxis(idip_,  hgrid%ax_modl, 2)
  call houaxis(idip_,  hgrid%ay_modl, 3)
  call houaxis(idip_,  hgrid%az_modl, 4)

  call from_param("nx_offd_pad_beg", hgrid%nx_offd_pad_beg, 50)
  call from_param("nx_offd_pad_end", hgrid%nx_offd_pad_end, 50)

  call from_param("ny_offd_pad_beg", hgrid%ny_offd_pad_beg, 0)
  call from_param("ny_offd_pad_end", hgrid%ny_offd_pad_end, 0)

  call from_param("nz_offd_pad_beg", hgrid%nz_offd_pad_beg, 50)
  call from_param("nz_offd_pad_end", hgrid%nz_offd_pad_end, 50) 

  call from_param("ntpx", ntpx, 0)
  call from_param("ntpy", ntpy, 0)
  call from_param("ntpz", ntpz, 0)

write(0,*)"hgrid%ax_offd%n,hgrid%ay_offd%n,hgrid%az_offd%n"
write(0,*) hgrid%ax_offd%n,hgrid%ay_offd%n,hgrid%az_offd%n
write(0,*)"hgrid%ax_modl%n,hgrid%ay_modl%n,hgrid%az_modl%n"
write(0,*) hgrid%ax_modl%n,hgrid%ay_modl%n,hgrid%az_modl%n
write(0,*)"hgrid%adip_agl%n", hgrid%adip_agl%n

  nkx = hgrid%ax_offd%n+hgrid%nx_offd_pad_beg+hgrid%nx_offd_pad_end
  nky = hgrid%ay_offd%n+hgrid%ny_offd_pad_beg+hgrid%ny_offd_pad_end
  nkz = hgrid%az_offd%n+hgrid%nz_offd_pad_beg+hgrid%nz_offd_pad_end

  allocate(hess(hgrid%ax_offd%n,hgrid%ay_offd%n,hgrid%az_offd%n))
  allocate(vmod(hgrid%ax_modl%n,hgrid%ay_modl%n,hgrid%az_modl%n))
  allocate(dip_agl(hgrid%adip_agl%n))
  allocate(sca_agl(hgrid%asca_agl%n))
  allocate(sca_map(nkx, nky, nkz))
  allocate(ahes(nkx, nky, nkz))

  call auxputch("n1", "i", nkx, test_)
  call auxputch("n2", "i", nky, test_)
  call auxputch("n3", "i", nkz, test_)
  call houaxis(test_, hgrid%ax_modl, 4)
  call houaxis(test_, hgrid%ay_modl, 5)
  call houaxis(test_, hgrid%az_modl, 6)

  call auxputch("n1", "i", nkx, msca_)
  call auxputch("n2", "i", nky, msca_)
  call auxputch("n3", "i", nkz, msca_)
  call houaxis(msca_, hgrid%ax_modl, 4)
  call houaxis(msca_, hgrid%ay_modl, 5)
  call houaxis(msca_, hgrid%az_modl, 6)


  hgrid%asca_agl%o = hgrid%asca_agl%o*PI/180. 
  hgrid%asca_agl%d = hgrid%asca_agl%d*PI/180.

  hgrid%adip_agl%o = hgrid%adip_agl%o*PI/180.
  hgrid%adip_agl%d = hgrid%adip_agl%d*PI/180.

  aphs%o = aphs%o*PI/180.
  aphs%d = aphs%d*PI/180.

  call sep_read(vmod, vmod_)

  stat = hes_gahter_init(hgrid, aphs, pct, ntpx, ntpy, ntpz, norm)

  do iz_modl=1, hgrid%az_modl%n; write(0,*)"iz_modl=", iz_modl
  do iy_modl=1, hgrid%ay_modl%n
  do ix_modl=1, hgrid%ax_modl%n
     vref = vmod(ix_modl,iy_modl,iz_modl)
     call hess_fileio_offd(hess, hgrid, ix_modl, iy_modl, iz_modl, hess_, "reed")
     stat = hes_gather_oper(hess, ahes, dip_agl, sca_agl, sca_map, vref, ws, ix_modl, iy_modl, iz_modl)
     call sep_write(ahes, test_)
     call sep_write(sca_map, msca_)  
     call sep_write(dip_agl, idip_)
     call sep_write(sca_agl, isca_)
  end do
  end do
  end do
end program
