module grn_fileio_mod
  use  base_struct_mod
  use  base_ioutil_mod
  use  sep
  implicit none

  integer, external, private :: sreed_window, srite_window, getch

contains
!================================================================
  subroutine grn_param_init(ggrid, pgrid, wblock, csou_, vmod_, wfds_)
    type(gre_grid)   :: ggrid
    type(pro_grid)   :: pgrid
    type(block)      :: wblock
    character(len=*) :: csou_, vmod_, wfds_

    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()

    csou_ = "csou"
    vmod_ = "vmod"
    wfds_ = "wfds"

    call from_param("radian", ggrid%radian, .true. )
    call from_param("spikes", ggrid%spikes, .false.)
    call from_param("ntpx"  , pgrid%ntpx  , 1      )
    call from_param("ntpy"  , pgrid%ntpy  , 1      )
    call from_param("nvrf"  , pgrid%nvrf  , 1      )
    call from_param("pdip"  , pgrid%pdip  , 65     )
    call from_param("oper"  , pgrid%oper  , "ssfpi")
    call from_param("trick" , pgrid%trick ,    0.14)
    call from_param("nws"   ,    wblock%ns, 1      )
    return
  end subroutine

  subroutine grn_histo_data(ggrid, pgrid, vgrid, csou_, vmod_)
    type(gre_grid)   :: ggrid
    type(pro_grid)   :: pgrid
    type(vel_grid)   :: vgrid
    character(len=*) :: csou_, vmod_

    real             :: surf_xmin_default, surf_xmax_default
    real             :: surf_ymin_default, surf_ymax_default 
    real             :: modl_xmin_default, modl_xmax_default
    real             :: modl_ymin_default, modl_ymax_default
    real             :: modl_zmin_default, modl_zmax_default

    if (.not.ggrid%spikes) then
       call auxin(csou_)
       call hinaxis(csou_, ggrid%aw, 1)
    else
       call from_param("ow", ggrid%aw%o, 0.)
       call from_param("dw", ggrid%aw%d, 1.)
       call from_param("nw", ggrid%aw%n, 1 )
    end if

    !!.. velocity parameter initialization
    call hinaxis(vmod_, vgrid%ax, 1)
    call hinaxis(vmod_, vgrid%ay, 2)
    call hinaxis(vmod_, vgrid%az, 3)

    !!.. input for acquisition parameters
    surf_xmin_default=vgrid%ax%o; surf_xmax_default=vgrid%ax%o+(vgrid%ax%n-1)*vgrid%ax%d
    surf_ymin_default=vgrid%ay%o; surf_ymax_default=vgrid%ay%o+(vgrid%ay%n-1)*vgrid%ay%d

    !!.. the receiver grids are assumed to be the same as the velocity grid
    call from_param("ox_surf", ggrid%ax_surf%o, vgrid%ax%o)
    call from_param("dx_surf", ggrid%ax_surf%d, vgrid%ax%d)
    call from_param("nx_surf", ggrid%ax_surf%n, vgrid%ax%n)

    call from_param("oy_surf", ggrid%ay_surf%o, vgrid%ay%o)
    call from_param("dy_surf", ggrid%ay_surf%d, vgrid%ay%d)
    call from_param("ny_surf", ggrid%ay_surf%n, vgrid%ay%n)

    ggrid%ax_surf%o=vgrid%ax%o+int(1.5+(ggrid%ax_surf%o-vgrid%ax%o)/vgrid%ax%d-1)*vgrid%ax%d
    ggrid%ay_surf%o=vgrid%ay%o+int(1.5+(ggrid%ay_surf%o-vgrid%ay%o)/vgrid%ay%d-1)*vgrid%ay%d

    ggrid%surf_xmin=ggrid%ax_surf%o; ggrid%surf_xmax=ggrid%ax_surf%o+(ggrid%ax_surf%n-1)*ggrid%ax_surf%d 
    ggrid%surf_ymin=ggrid%ay_surf%o; ggrid%surf_ymax=ggrid%ay_surf%o+(ggrid%ay_surf%n-1)*ggrid%ay_surf%d

    !!.. input for model points, the samplings are assumed to be the same as
    !!.. those of the velocity grid

    modl_xmin_default=vgrid%ax%o; modl_xmax_default=vgrid%ax%o+(vgrid%ax%n-1)*vgrid%ax%d
    modl_ymin_default=vgrid%ay%o; modl_ymax_default=vgrid%ay%o+(vgrid%ay%n-1)*vgrid%ay%d
    modl_zmin_default=vgrid%az%o; modl_zmax_default=vgrid%az%o+(vgrid%az%n-1)*vgrid%az%d

    call from_param("modl_xmin", ggrid%modl_xmin, modl_xmin_default); call from_param("modl_xmax", ggrid%modl_xmax, modl_xmax_default)
    call from_param("modl_ymin", ggrid%modl_ymin, modl_ymin_default); call from_param("modl_ymax", ggrid%modl_ymax, modl_ymax_default)
    call from_param("modl_zmin", ggrid%modl_zmin, modl_zmin_default); call from_param("modl_zmax", ggrid%modl_zmax, modl_zmax_default)

    ggrid%modl_xmin=vgrid%ax%o+int(1.5+(ggrid%modl_xmin-vgrid%ax%o)/vgrid%ax%d-1)*vgrid%ax%d
    ggrid%modl_ymin=vgrid%ay%o+int(1.5+(ggrid%modl_ymin-vgrid%ay%o)/vgrid%ay%d-1)*vgrid%ay%d
    ggrid%modl_zmin=vgrid%az%o+int(1.5+(ggrid%modl_zmin-vgrid%az%o)/vgrid%az%d-1)*vgrid%az%d

    ggrid%modl_xmax=vgrid%ax%o+int(1.5+(ggrid%modl_xmax-vgrid%ax%o)/vgrid%ax%d-1)*vgrid%ax%d
    ggrid%modl_ymax=vgrid%ay%o+int(1.5+(ggrid%modl_ymax-vgrid%ay%o)/vgrid%ay%d-1)*vgrid%ay%d
    ggrid%modl_zmax=vgrid%az%o+int(1.5+(ggrid%modl_zmax-vgrid%az%o)/vgrid%az%d-1)*vgrid%az%d

    ggrid%ax_modl%o=ggrid%modl_xmin; ggrid%ax_modl%d=vgrid%ax%d; ggrid%ax_modl%n=1.5+(ggrid%modl_xmax-ggrid%modl_xmin)/ggrid%ax_modl%d
    ggrid%ay_modl%o=ggrid%modl_ymin; ggrid%ay_modl%d=vgrid%ay%d; ggrid%ay_modl%n=1.5+(ggrid%modl_ymax-ggrid%modl_ymin)/ggrid%ay_modl%d
    ggrid%az_modl%o=ggrid%modl_zmin; ggrid%az_modl%d=vgrid%az%d; ggrid%az_modl%n=1.5+(ggrid%modl_zmax-ggrid%modl_zmin)/ggrid%az_modl%d

    !!.. input padding
    call from_param("nxpad_beg", pgrid%nxpad_beg, 0); call from_param("nxpad_end", pgrid%nxpad_end, 0)
    call from_param("nypad_beg", pgrid%nypad_beg, 0); call from_param("nypad_end", pgrid%nypad_end, 0)

    return 
  end subroutine

  subroutine grn_histo_wfds(ggrid, wfds_)
    type(gre_grid)   :: ggrid
    character(len=*) :: wfds_

    call auxout(wfds_)

    call houaxis(wfds_, ggrid%ax_modl, 1)
    call houaxis(wfds_, ggrid%ay_modl, 2)
    call houaxis(wfds_, ggrid%az_modl, 3)
    call houaxis(wfds_, ggrid%aw     , 4)
    call houaxis(wfds_, ggrid%ax_surf, 5)
    call houaxis(wfds_, ggrid%ay_surf, 6)

    call auxputch("esize" , "i",         8, wfds_)
    call auxputch("label1", "s", "x_modl" , wfds_)
    call auxputch("label2", "s", "y_modl" , wfds_)
    call auxputch("label3", "s", "z_modl" , wfds_)
    call auxputch("label4", "s", "freq."  , wfds_)
    call auxputch("label5", "s", "x_surf" , wfds_)
    call auxputch("label6", "s", "y_surf" , wfds_)

    return
  end subroutine

!=============================================================================================
  subroutine grn_fio_csou(csou, ggrid, indw, nws, filetag)
    complex, dimension(:) :: csou
    type(gre_grid)        :: ggrid
    integer               :: indw, nws
    character(len=*)      :: filetag

    integer, dimension(1) :: n_cube, j_wind, n_wind, f_wind
    integer               :: n_dim_cube, err

    n_dim_cube = 1
    n_cube = ggrid%aw%n
    j_wind = 1
    n_wind = nws
    f_wind = indw - 1
    csou = 0.
    err = sreed_window(filetag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,csou)
    if (err.lt.0) call seperr("wei_greens_io_csou: trouble reading source wavefield")
  end subroutine
!=====================================================================================
  subroutine grn_fio_velo(vmod, vgrid, pgrid, indz, nzs, filetag)
    real, dimension(:,:,:) :: vmod
    type(vel_grid)         :: vgrid
    type(pro_grid)         :: pgrid
    integer                :: indz
    integer                :: nzs
    character(len=*)       :: filetag
    integer, dimension(3)  :: n_cube, j_wind, n_wind, f_wind
    integer                :: n_dim_cube, err

    n_dim_cube = 3
    n_cube = (/vgrid%ax%n, vgrid%ay%n, vgrid%az%n /)
    j_wind = (/1, 1, 1/)
    n_wind = (/pgrid%ax%n, pgrid%ay%n, nzs/)
    f_wind = (/vgrid%index_provel_xmin-1, vgrid%index_provel_ymin-1, indz-1/)

    vmod = 0.
    err = sreed_window(filetag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,vmod)
    if (err .lt. 0) call seperr("sreed_window: trouble reading velocity file")
    return
  end subroutine

  subroutine grn_fio_wfds_init(ggrid, filetag)
    type(gre_grid)                           :: ggrid
    character(len=*)                         :: filetag

    complex, dimension(:,:,:,:), allocatable :: tmp
    integer                                  :: ix_surf, iy_surf

    allocate(tmp(ggrid%ax_modl%n,ggrid%ay_modl%n,ggrid%az_modl%n,ggrid%aw%n))
    tmp = 0.
    do iy_surf=1, ggrid%ay_surf%n
    do ix_surf=1, ggrid%ax_surf%n
       call sep_write(tmp, filetag)
    end do
    end do
    return
  end subroutine
!================================================================================
! write
!--------------------------------------------------------------------------------  
  subroutine grn_fio_wfds(wfds, ggrid, ind_ws, nws, ind_surfx, ind_surfy, filetag, flag)
    complex, dimension(:,:,:,:), target  :: wfds
    type(gre_grid)                       :: ggrid
    integer                              :: nws, ind_ws, ind_surfx, ind_surfy
    character(len=*)                     :: filetag
    character(len=4)                     :: flag

    integer, dimension(6)                :: n_cube, j_wind, n_wind, f_wind
    integer                              :: n_dim_cube, err
    integer                              :: indx_xmin, indx_xmax
    integer                              :: indx_ymin, indx_ymax
    integer                              :: indx_zmin, indx_zmax
    complex, dimension(:,:,:,:), pointer :: pwfds

    n_dim_cube = 6
    n_cube = (/ggrid%ax_modl%n,ggrid%ay_modl%n,ggrid%az_modl%n,ggrid%aw%n,ggrid%ax_surf%n,ggrid%ay_surf%n/)
    j_wind = (/1, 1, 1, 1, 1, 1/)
    n_wind = (/ggrid%ax_modl%n,ggrid%ay_modl%n,ggrid%az_modl%n,nws, 1, 1/)
    f_wind = (/0, 0, 0, ind_ws-1, ind_surfx-1, ind_surfy-1/)

    indx_xmin = ggrid%index_modpro_xmin; indx_xmax = ggrid%index_modpro_xmax
    indx_ymin = ggrid%index_modpro_ymin; indx_ymax = ggrid%index_modpro_ymax
    indx_zmin = ggrid%index_modpro_zmin; indx_zmax = ggrid%index_modpro_zmax

    pwfds => wfds(indx_xmin:indx_xmax,indx_ymin:indx_ymax,indx_zmin:indx_zmax,:)

    if (flag == "rite") then
       err = srite_window(filetag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,pwfds)
    else if(flag == "reed") then
       pwfds = 0.
       err = sreed_window(filetag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,pwfds)
    else
       call seperr("you have to specify the flag for wei_io_greens()")
    end if
    nullify(pwfds)
    return
  end subroutine


end module








