module wem_fileio_mod
  use  base_struct_mod
  use  base_ioutil_mod
  use  sep
  implicit none

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

contains
!================================================================
  subroutine wem_param_init(dgrid, pgrid, igrid, hgrid, wblock, crec_, csou_, vmod_, imag_, hess_, hsou_)
    type(dat_grid)   :: dgrid
    type(pro_grid)   :: pgrid
    type(img_grid)   :: igrid
    type(hes_grid)   :: hgrid
    type(block)      :: wblock
    character(len=*) :: crec_, csou_, vmod_, imag_, hess_, hsou_

    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" 
    crec_ = "crec" 
    vmod_ = "vmod" 
    imag_ = "imag"
    hess_ = "hess"
    hsou_ = "hsou"

    call from_param("radian", dgrid%radian, .true. )
    call from_param("offset", dgrid%offset, .false.)
    call from_param("comsou", dgrid%comsou, .false.)
    call from_param("hdrsou", dgrid%hdrsou, .false.)
    call from_param("wntimg", igrid%wntimg, .true. )
    call from_param("wnthes", hgrid%wnthes, .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("hestyp", hgrid%hestyp, 1      )
    call from_param("nws"   ,    wblock%ns, 1      )
    return
  end subroutine
!================================================================
! read parameters from histories
!----------------------------------------------------------------    
  subroutine wem_histo_data(dgrid, vgrid, igrid, pgrid, crec_, csou_, vmod_, hsou_)
    type(dat_grid)      :: dgrid
    type(vel_grid)      :: vgrid
    type(img_grid)      :: igrid
    type(pro_grid)      :: pgrid
    character(len=*)    :: crec_, csou_, vmod_, hsou_
    type(axis)          :: agxcs, agycs, asxcs, asycs 
    real                :: image_xmin_default, image_xmax_default, image_dx_default
    real                :: image_ymin_default, image_ymax_default, image_dy_default
    real                :: image_zmin_default, image_zmax_default, image_dz_default
    real                :: ohx_default, ohy_default
    integer             :: nshg

    call auxin(crec_)
    call auxin(csou_)
    call auxin(vmod_)

    !!.. read the receiver coordinates, which are assumed to be regular
    call hinaxis(crec_, dgrid%agx, 1)
    call hinaxis(crec_, dgrid%agy, 2)
    call hinaxis(crec_, dgrid%aw , 3)
    !!.. total receiver traces for one shot
    dgrid%ntrg = dgrid%agx%n * dgrid%agy%n
    !!.. if source header is presented
    if (dgrid%hdrsou) then
       call auxin(hsou_)
       call from_aux(hsou_, "n2", dgrid%nshg)   !!.. total shots
       call from_aux(crec_, "n4", nshg)
       if (dgrid%nshg /= nshg) call seperr("ERR: hsou-nshg /= crec-nshg")
    else
       call hinaxis(crec_, dgrid%asx, 4)
       call hinaxis(crec_, dgrid%asy, 5)
       dgrid%nshg = dgrid%asx%n * dgrid%asy%n   !!.. total shots
    end if
    !!.. if composite sources are used 
    if (dgrid%comsou) then
       call hinaxis(csou_, agxcs, 1)
       call hinaxis(csou_, agycs, 2)
       call hinaxis(csou_, asxcs, 4)
       call hinaxis(csou_, asycs, 5)
       dgrid%ntrs = agxcs%n * agycs%n           !!.. size of composite source
       dgrid%nshs = asxcs%n * asycs%n           !!.. number of composite sources
       if (dgrid%ntrs /= dgrid%ntrg) call seperr("ERROR: dgrid%ntrs /= dgrid%ntrg")
       if (dgrid%nshs /= dgrid%nshg) call seperr("ERROR: dgrid%nshs /= dgrid%nshg")
    else
       dgrid%ntrs = 1
       dgrid%nshs = 1
    end if

    if (dgrid%offset) then
       dgrid%gxmin_save = dgrid%agx%o
       dgrid%gymin_save = dgrid%agy%o
    end if

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

    !!.. the default image size is set to be the same as the velocity model
    image_xmin_default = vgrid%ax%o
    image_xmax_default = vgrid%ax%o + (vgrid%ax%n-1)*vgrid%ax%d
    image_dx_default   = vgrid%ax%d

    image_ymin_default = vgrid%ay%o
    image_ymax_default = vgrid%ay%o + (vgrid%ay%n-1)*vgrid%ay%d
    image_dy_default   = vgrid%ay%d

    image_zmin_default = vgrid%az%o
    image_zmax_default = vgrid%az%o + (vgrid%az%n-1)*vgrid%az%d
    image_dz_default   = vgrid%az%d

    call from_param("image_xmin", igrid%image_xmin, image_xmin_default)
    call from_param("image_xmax", igrid%image_xmax, image_xmax_default)
    call from_param("image_ymin", igrid%image_ymin, image_ymin_default)
    call from_param("image_ymax", igrid%image_ymax, image_ymax_default)
    call from_param("image_zmin", igrid%image_zmin, image_zmin_default)
    call from_param("image_zmax", igrid%image_zmax, image_zmax_default)
    call from_param("image_dx"  , igrid%image_dx  , image_dx_default  )
    call from_param("image_dy"  , igrid%image_dy  , image_dy_default  )
    call from_param("image_dz"  , igrid%image_dz  , image_dz_default  )

    igrid%image_xmin = vgrid%ax%o + int(1.5+(igrid%image_xmin-vgrid%ax%o)/vgrid%ax%d-1)*vgrid%ax%d
    igrid%image_ymin = vgrid%ay%o + int(1.5+(igrid%image_ymin-vgrid%ay%o)/vgrid%ay%d-1)*vgrid%ay%d
    igrid%image_zmin = vgrid%az%o + int(1.5+(igrid%image_zmin-vgrid%az%o)/vgrid%az%d-1)*vgrid%az%d

    igrid%image_xmax = vgrid%ax%o + int(1.5+(igrid%image_xmax-vgrid%ax%o)/vgrid%ax%d-1)*vgrid%ax%d
    igrid%image_ymax = vgrid%ay%o + int(1.5+(igrid%image_ymax-vgrid%ay%o)/vgrid%ay%d-1)*vgrid%ay%d
    igrid%image_zmax = vgrid%az%o + int(1.5+(igrid%image_zmax-vgrid%az%o)/vgrid%az%d-1)*vgrid%az%d

    if(igrid%image_xmin < image_xmin_default) call seperr("image_xmin is out of the input velocity model")
    if(igrid%image_xmax > image_xmax_default) call seperr("image_xmax is out of the input velocity model")
    if(igrid%image_ymin < image_ymin_default) call seperr("image_ymin is out of the input velocity model")
    if(igrid%image_ymax > image_ymax_default) call seperr("image_ymax is out of the input velocity model")
    if(igrid%image_zmin < image_zmin_default) call seperr("image_zmin is out of the input velocity model")
    if(igrid%image_zmax > image_zmax_default) call seperr("image_zmax is out of the input velocity model")

    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)
    call from_param("nzpad_beg", pgrid%nzpad_beg, 0)
    call from_param("nzpad_end", pgrid%nzpad_end, 0)

    call from_param("nhx", igrid%ahx%n, 1         )
    call from_param("nhy", igrid%ahy%n, 1         )
    call from_param("nhz", igrid%ahz%n, 1         )

    call from_param("dhx", igrid%ahx%d, vgrid%ax%d)
    call from_param("dhy", igrid%ahy%d, vgrid%ay%d)
    call from_param("dhz", igrid%ahz%d, vgrid%az%d)

    call from_param("ohx", igrid%ahx%o, 0.        ) 
    call from_param("ohy", igrid%ahy%o, 0.        )
    call from_param("ohz", igrid%ahz%o, 0.        )
    return
  end subroutine

!===========================================================================
! create image history
!  adj=F, input history file for reading data
!  adj=T, output history file for writing data
!---------------------------------------------------------------------------
  subroutine wem_histo_imag(igrid, imag_)
    type(img_grid)   :: igrid
    character(len=*) :: imag_

    !!.. define the image grid
    igrid%ax%o = igrid%image_xmin
    igrid%ax%d = igrid%image_dx
    igrid%ax%n = 1.5 + (igrid%image_xmax - igrid%image_xmin)/igrid%image_dx

    igrid%ay%o = igrid%image_ymin
    igrid%ay%d = igrid%image_dy
    igrid%ay%n = 1.5 + (igrid%image_ymax - igrid%image_ymin)/igrid%image_dy

    igrid%az%o = igrid%image_zmin
    igrid%az%d = igrid%image_dz
    igrid%az%n = 1.5 + (igrid%image_zmax - igrid%image_zmin)/igrid%image_dz 

    call auxinout(imag_)

    call houaxis(imag_, igrid%ax , 1)
    call houaxis(imag_, igrid%ay , 2)
    call houaxis(imag_, igrid%az , 3)
    call houaxis(imag_, igrid%ahx, 4)
    call houaxis(imag_, igrid%ahy, 5)
    call houaxis(imag_, igrid%ahz, 6)

    call auxputch("label1", "s", "x" , imag_)
    call auxputch("label2", "s", "y" , imag_)
    call auxputch("label3", "s", "z" , imag_)
    call auxputch("label4", "s", "hx", imag_)
    call auxputch("label5", "s", "hy", imag_)
    call auxputch("label6", "s", "hz", imag_)

    return
  end subroutine 
!--------------------------------------------------------------------------------------
  subroutine wem_histo_hess(hgrid, igrid, hess_)
    type(hes_grid)   :: hgrid
    type(img_grid)   :: igrid
    character(len=*) :: hess_
    hgrid%ax_modl = igrid%ax
    hgrid%ay_modl = igrid%ay
    hgrid%az_modl = igrid%az

    call from_param("nx_offd", hgrid%ax_offd%n, 1 )
    call from_param("ox_offd", hgrid%ax_offd%o, 0.)
    call from_param("dx_offd", hgrid%ax_offd%d, 1.)

    call from_param("ny_offd", hgrid%ay_offd%n, 1 )
    call from_param("oy_offd", hgrid%ay_offd%o, 0.)
    call from_param("dy_offd", hgrid%ay_offd%d, 1.)

    call from_param("nz_offd", hgrid%az_offd%n, 1 )
    call from_param("oz_offd", hgrid%az_offd%o, 0.)
    call from_param("dz_offd", hgrid%az_offd%d, 1.)
 
    call auxinout(hess_)

    call houaxis(hess_, hgrid%ax_modl, 1)
    call houaxis(hess_, hgrid%ay_modl, 2)
    call houaxis(hess_, hgrid%az_modl, 3)
    call houaxis(hess_, hgrid%ax_offd, 4) 
    call houaxis(hess_, hgrid%ay_offd, 5)
    call houaxis(hess_, hgrid%az_offd, 6)
    call auxputch("label1", "s", "x_modl", hess_)
    call auxputch("label2", "s", "y_modl", hess_)
    call auxputch("label3", "s", "z_modl", hess_)
    call auxputch("label4", "s", "x_offd", hess_)
    call auxputch("label5", "s", "y_offd", hess_)
    call auxputch("label6", "s", "z_offd", hess_)

    return
  end subroutine
!--------------------------------------------------------------------------------------
  subroutine wem_fio_hess_init(filetag, hgrid)
    character(len=*)                :: filetag
    type(hes_grid)                  :: hgrid
    real, allocatable, dimension(:) :: tmp
    integer                         :: ioffd, nmodl, noffd
    nmodl = hgrid%ax_modl%n*hgrid%ay_modl%n*hgrid%az_modl%n
    noffd = hgrid%ax_offd%n*hgrid%ay_offd%n*hgrid%az_offd%n
    allocate(tmp(nmodl))
    tmp = 0.
    do ioffd=1, noffd
       call sep_write(tmp, filetag, esize=4)
    end do
    deallocate(tmp)
    return
  end subroutine
!-------------------------------------------------------------------------------------
  subroutine wem_fio_imag_init(filetag, igrid)
    character(len=*)                :: filetag
    type(img_grid)                  :: igrid

    real, allocatable, dimension(:) :: tmp
    integer                         :: ioffs, nmodl, noffs
    nmodl = igrid%ax%n*igrid%ay%n*igrid%az%n
    noffs = igrid%ahx%n*igrid%ahy%n*igrid%ahz%n
    allocate(tmp(nmodl))
    tmp = 0.
    do ioffs=1, noffs
       call sep_write(tmp, filetag, esize=4)
    end do
    deallocate(tmp)
    return
  end subroutine

!=====================================================================================
! read or write one depth block image
!-------------------------------------------------------------------------------------
   subroutine wem_fio_imag(modl, igrid, ogrid, indz, nzs, filetag, flag)
     real, dimension(:,:,:,:,:) :: modl
     type(img_grid)             :: igrid
     type(out_grid)             :: ogrid
     integer                    :: indz
     integer                    :: nzs
     character(len=*)                    :: filetag
     character(len=4)                    :: flag
            
     integer, dimension(5)               :: n_cube, j_wind, n_wind, f_wind
     integer                             :: n_dim_cube, err


     n_dim_cube = 5
     n_cube = (/igrid%ax%n, igrid%ay%n, igrid%ahx%n, igrid%ahy%n, igrid%az%n /)
     j_wind = (/1, 1, 1, 1, 1/)
     n_wind = (/ogrid%ax%n, ogrid%ay%n, ogrid%ahx%n, ogrid%ahy%n, nzs/)
     f_wind = (/ogrid%index_outimg_xmin-1, ogrid%index_outimg_ymin-1, 0, 0, indz-1/)

     if (flag == "rite") then
        err = srite_window(filetag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,modl)
        if (err .lt. 0) call seperr("wem_io_imag: trouble writing image file")
     else if(flag == "reed") then
        modl = 0.
        err = sreed_window(filetag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,modl)
        if (err .lt. 0) call seperr("wem_io_imag: trouble reading image file")
     else
        call seperr("you have to specify the flag for wem_io_imag()")
     end if
     return
   end subroutine

!=============================================================================
! read one frequency block receiver wavefield
!-----------------------------------------------------------------------------
   subroutine wem_fio_crec(crec, dgrid, indw, nws, ish, filetag, flag)
     complex, dimension(:,:) :: crec
     type(dat_grid)          :: dgrid 
     integer                 :: indw, nws, ish
     character(len=*)        :: filetag
     character(len=4)        :: flag
     integer, dimension(3)   :: n_cube, j_wind, n_wind, f_wind
     integer                 :: n_dim_cube, err  

     n_dim_cube = 3
     n_cube = (/dgrid%ntrg, dgrid%aw%n, dgrid%nshg/)
     j_wind = (/1, 1, 1/)
     n_wind = (/dgrid%ntrg, nws, 1/)
     f_wind = (/0, indw-1, ish-1/)

     if (flag == "rite") then
        err = srite_window(filetag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,crec)
        if (err.lt.0) call seperr("wem_fio_crec: trouble writing receiver wavefield")
     else if(flag == "reed") then
        crec = 0.
        err = sreed_window(filetag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,crec)
        if (err.lt.0) call seperr("wem_fio_crec: trouble reading receiver wavefield")
     else
        call seperr("you have to specify the flag for wem_fio_crec()")
     end if
     return
   end subroutine

!=============================================================================
! read one frequency block receiver wavefield
!-----------------------------------------------------------------------------
   subroutine wem_fio_csou(csou, dgrid, indw, nws, ish, filetag, flag)
     complex, dimension(:,:) :: csou
     type(dat_grid)          :: dgrid
     integer                 :: indw, nws, ish
     character(len=*)        :: filetag
     character(len=4)        :: flag
     integer, dimension(3)   :: n_cube, j_wind, n_wind, f_wind
     integer                 :: n_dim_cube, err

     n_dim_cube = 3
     n_cube = (/dgrid%ntrs, dgrid%aw%n, dgrid%nshs/)
     j_wind = (/1, 1, 1/)
     n_wind = (/dgrid%ntrs, nws, 1/)
     if (dgrid%comsou) then 
        f_wind = (/0, indw-1, ish-1/)
     else
        f_wind = (/0, indw-1, 0/)
     end if   
     if (flag == "rite") then
        err = srite_window(filetag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,csou)
        if (err.lt.0) call seperr("wem_fio_csou: trouble writing receiver wavefield")
     else if(flag == "reed") then
        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("wem_fio_csou: trouble reading receiver wavefield")
     else
        call seperr("you have to specify the flag for wem_fio_csou()")
     end if
     return
   end subroutine

!==============================================================================
! read one depth block velocity model
!------------------------------------------------------------------------------
  subroutine wem_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 wem_wfds_alloc(nws)
!    integer :: nws
!    if (dgrid%comsou) then
!       allocate(csou(dgrid%ntrs, nws))
!    else 
!       allocate(csou(1, nws))
!    end if
!    allocate(wfds(pgrid%ax%n, pgrid%ay%n, pgrid%nzstep, nws))
!    wfds = cmplx(0., 0.)
!    csou = cmplx(0., 0.)
!    return
!  end subroutine
!  subroutine wem_wfds_dealloc()
!    deallocate(csou)
!    deallocate(wfds)
!    return
!  end subroutine
!
!  subroutine wem_wfdr_alloc(nws)
!    integer :: nws
!    allocate(crec(dgrid%ntrg, nws))
!    allocate(wfdr(pgrid%ax%n, pgrid%ay%n, pgrid%nzstep, nws))
!    crec = cmplx(0., 0.)
!    wfdr = cmplx(0., 0.)
!    return
!  end subroutine
!  subroutine wem_wfdr_dealloc()
!    deallocate(crec)
!    deallocate(wfdr)
!    return
!  end subroutine
!
!  subroutine wem_modl_alloc()
!    integer :: nzs
!    allocate(modl(ogrid%ax%n, ogrid%ay%n, ogrid%az%n))
!    modl = 0.
!    return
!  end subroutine
!  subroutine wem_modl_dealloc()
!    deallocate(modl)
!    return
!  end subroutine
!
!  subroutine wem_vmod_alloc(nzs)
!    integer :: nzs
!    allocate(vmod(pgrid%ax%n, pgrid%ay%n, nzs))
!    allocate(vref(nvrf, nzs))
!    vmod = 0.
!    vref = 0.
!    return
!  end subroutine
!  subroutine wem_vmod_dealloc()
!    deallocate(vmod)
!    deallocate(vref)
!    return
!  end subroutine
!
!  subroutine wem_hess_alloc()
!    allocate(hess(ogrid%ax%n, ogrid%ay%n, ogrid%az%n))
!    hess = 0.
!    return
!  end subroutine
!  subroutine wem_hess_dealloc()
!    deallocate(hess)
!    return
!  end subroutine
end module














