module wei_io_mod
  use  wei_global_defi_mod
  use  sep
  implicit none

  integer, external :: sreed_window, srite_window, getch

contains
!================================================================
  subroutine wei_param_init()
    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()

    sou_     = "sou"    ; rec_     = "rec"   ; vel_    = "vel"   ; img_ = "img"
    wfld_d_  = "wfld_d" ; wfld_u_  = "wfld_u"
    hdr_sou_ = "hdr_sou"; hdr_rec_ = "hdr_rec"
    hessian_ = "hessian"; greens_  = "greens"; datumn_ = "datumn"

    call from_param("adj"         , adj          , .true. )
    call from_param("offset"      , offset       , .false.)
    call from_param("nref"        , nref         , 1      )
    call from_param("lloyd"       , lloyd        , .false.)
    call from_param("decon"       , decon        , .false.)
    call from_param("radian"      , radian       , .false.)
    call from_param("taper"       , taper        , .false.)
    call from_param("taper_pct"   , taper_pct    , 0.1    )
    call from_param("taper_type"  , taper_type   , 1      )
    call from_param("num_report"  , num_report   , 5      )
    call from_param("eps"         , eps          , 0.001  )
    call from_param("maxang"      , maxang       , 90.    )
    call from_param("hessian_type", hessian_type , 0      )
    call from_param("greens"      , greens       , .false.)
    call from_param("wfld_d_out"  , wfld_d_out   , .false.)
    call from_param("wfld_u_out"  , wfld_u_out   , .false.)
  end subroutine
!================================================================
! read parameters from histories
!----------------------------------------------------------------    
  subroutine wei_data_history(adj)
    logical, intent(in) :: adj
    real                :: image_xmin_default, image_xmax_default
    real                :: image_ymin_default, image_ymax_default
    real                :: image_dx_default  , image_dy_default
    real                :: ohx_default, ohy_default
    integer             :: num_match
    character(len=1024) :: string


    dgrid%hdr_rec = .false.
    dgrid%hdr_sou = .false.

    if (adj) then !! in case of migration

       call auxin(rec_)
       call auxin(sou_)
       call auxin(vel_)
       call auxinout(wfld_d_)
       if(wfld_u_out) call auxout(wfld_u_)
       !!.. if providing a receiver header, a source header should also be provided
       num_match = getch(hdr_rec_, "s", string) + getch(hdr_sou_, "s", string)
       if (num_match == 2) dgrid%hdr_rec = .true.

       num_match = getch(hdr_sou_, "s", string)
       if (num_match == 1) dgrid%hdr_sou = .true.

       !!.. if the receiver header is presented
       if (dgrid%hdr_rec) then
          call from_aux(rec_, "n1", dgrid%ntr)
          call from_aux(rec_, "n2", dgrid%aw%n)
          call from_aux(rec_, "d2", dgrid%aw%d)
          call from_aux(rec_, "o2", dgrid%aw%o)
          call auxin(hdr_rec_)
       !!.. if the receiver header is not presented
       else
          call from_aux(rec_, "n1", dgrid%agx%n)
          call from_aux(rec_, "d1", dgrid%agx%d)
          call from_aux(rec_, "o1", dgrid%agx%o)

          call from_aux(rec_, "n2", dgrid%agy%n)
          call from_aux(rec_, "d2", dgrid%agy%d)
          call from_aux(rec_, "o2", dgrid%agy%o)

          call from_aux(rec_, "n3", dgrid%aw%n)
          call from_aux(rec_, "d3", dgrid%aw%d)
          call from_aux(rec_, "o3", dgrid%aw%o)
          dgrid%ntr = dgrid%agx%n * dgrid%agy%n
       end if

       !!.. if the source header is presented
       if (dgrid%hdr_sou) then
          call from_aux(rec_, "n4", dgrid%nsh)
          call auxin(hdr_sou_)
       !!.. if the source header is not presented
       else
          call from_aux(rec_, "n4", dgrid%asx%n)
          call from_aux(rec_, "d4", dgrid%asx%d)
          call from_aux(rec_, "o4", dgrid%asx%o)

          call from_aux(rec_, "n5", dgrid%asy%n)
          call from_aux(rec_, "d5", dgrid%asy%d)
          call from_aux(rec_, "o5", dgrid%asy%o)
          dgrid%nsh = dgrid%asx%n * dgrid%asy%n
       end if

       !!.. velocity parameter initialization
       call from_aux(vel_, "n1", vgrid%ax%n)
       call from_aux(vel_, "d1", vgrid%ax%d)
       call from_aux(vel_, "o1", vgrid%ax%o)

       call from_aux(vel_, "n2", vgrid%ay%n)
       call from_aux(vel_, "d2", vgrid%ay%d)
       call from_aux(vel_, "o2", vgrid%ay%o)

       call from_aux(vel_, "n3", vgrid%az%n)
       call from_aux(vel_, "d3", vgrid%az%d)
       call from_aux(vel_, "o3", vgrid%az%o)

       !vgrid%az%n = vgrid%az%n - 1
       !vgrid%az%o = vgrid%az%o + vgrid%az%d

       !!.. read parameters for freq and depth blocks
       call from_param("nws", wblock%ns, dgrid%aw%n)
       call from_param("nzs", zblock%ns, vgrid%az%n)

       !!.. 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

       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_dx"  , igrid%image_dx  , image_dx_default)
       call from_param("image_dy"  , igrid%image_dy  , image_dy_default)

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

       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("jpx", pgrid%jx, 1.)
       call from_param("jpy", pgrid%jy, 1.)

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

       if (mod(igrid%ahx%n, 2) == 0) igrid%ahx%n = igrid%ahx%n + 1
       if (mod(igrid%ahy%n, 2) == 0) igrid%ahy%n = igrid%ahy%n + 1

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

       ohx_default = - int(igrid%ahx%n/2) * igrid%ahx%d
       ohy_default = - int(igrid%ahy%n/2) * igrid%ahy%d

       call from_param("ohx", igrid%ahx%o, ohx_default) 
       call from_param("ohy", igrid%ahy%o, ohy_default)


    else !! in case of demigration

       call auxin(sou_)
       call auxin(vel_)
       call auxinout(rec_)
       call auxinout(wfld_d_)
       !!.. if providing a receiver header, a source header should also be provided
       num_match = getch(hdr_rec_, "s", string) + getch(hdr_sou_, "s", string)
       if (num_match == 2) dgrid%hdr_rec = .true.

       num_match = getch(hdr_sou_, "s", string)
       if (num_match == 1) dgrid%hdr_sou = .true.

       call from_aux(sou_, "n3", dgrid%aw%n)
       call from_aux(sou_, "d3", dgrid%aw%d)
       call from_aux(sou_, "o3", dgrid%aw%o)

       !!.. velocity parameter initialization
       call from_aux(vel_, "n1", vgrid%ax%n)
       call from_aux(vel_, "d1", vgrid%ax%d)
       call from_aux(vel_, "o1", vgrid%ax%o)

       call from_aux(vel_, "n2", vgrid%ay%n)
       call from_aux(vel_, "d2", vgrid%ay%d)
       call from_aux(vel_, "o2", vgrid%ay%o)

       call from_aux(vel_, "n3", vgrid%az%n)
       call from_aux(vel_, "d3", vgrid%az%d)
       call from_aux(vel_, "o3", vgrid%az%o)

       !vgrid%az%n = vgrid%az%n - 1
       !vgrid%az%o = vgrid%az%o + vgrid%az%d

       !!.. read parameters for freq and depth blocks
       call from_param("nws", wblock%ns, dgrid%aw%n)
       call from_param("nzs", zblock%ns, vgrid%az%n)

       !!.. 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

       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_dx"  , igrid%image_dx  , image_dx_default)
       call from_param("image_dy"  , igrid%image_dy  , image_dy_default)

       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 the receiver header is provided for receiver geometry
       if (dgrid%hdr_rec) then
          call auxin(hdr_rec_)
          call from_aux(hdr_rec_, "n2", dgrid%ntr)
          call seperr("not working yet")
       else !!.. if the receiver header is not provided
          call from_param("ngx", dgrid%agx%n, 1)
          call from_param("dgx", dgrid%agx%d, 1.)
          call from_param("ogx", dgrid%agx%o, 0.)

          call from_param("ngy", dgrid%agy%n, 1)
          call from_param("dgy", dgrid%agy%d, 1.)
          call from_param("ogy", dgrid%agy%o, 0.)
          dgrid%ntr = dgrid%agx%n * dgrid%agy%n
       end if

       !!.. if the source header is presented
       if (dgrid%hdr_sou) then
          call auxin(hdr_sou_)
          call from_aux(hdr_sou_, "n2", dgrid%nsh)
       else !!.. if the source header is not presented
          call from_param("nsx", dgrid%asx%n, 1)
          call from_param("dsx", dgrid%asx%d, 1.)
          call from_param("osx", dgrid%asx%o, 0.)
     
          call from_param("nsy", dgrid%asy%n, 1)
          call from_param("dsy", dgrid%asy%d, 1.)
          call from_param("osy", dgrid%asy%o, 0.)
          dgrid%nsh = dgrid%asx%n * dgrid%asy%n
       end if

       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("jpx", pgrid%jx, 1.)
       call from_param("jpy", pgrid%jy, 1.)

       !!.. the data header
       call auxputch("n1", "i", dgrid%agx%n, rec_)
       call auxputch("d1", "r", dgrid%agx%d, rec_)
       call auxputch("o1", "r", dgrid%agx%o, rec_)

       call auxputch("n2", "i", dgrid%agy%n, rec_)
       call auxputch("d2", "r", dgrid%agy%d, rec_)
       call auxputch("o2", "r", dgrid%agy%o, rec_)

       call auxputch("n3", "i", dgrid%aw%n, rec_)
       call auxputch("d3", "r", dgrid%aw%d, rec_)
       call auxputch("o3", "r", dgrid%aw%o, rec_)

       call auxputch("n4", "i", dgrid%asx%n, rec_)
       call auxputch("d4", "r", dgrid%asx%d, rec_)
       call auxputch("o4", "r", dgrid%asx%o, rec_)

       call auxputch("n5", "i", dgrid%asy%n, rec_)
       call auxputch("d5", "r", dgrid%asy%d, rec_)
       call auxputch("o5", "r", dgrid%asy%o, rec_)

       call auxputch("esize" , "i",  8    , rec_)
       call auxputch("label1", "s", "gx"  , rec_)
       call auxputch("label2", "s", "gy"  , rec_)
       call auxputch("label3", "s", "freq", rec_)
       call auxputch("label4", "s", "sx"  , rec_)
       call auxputch("label5", "s", "sy"  , rec_)
    end if
  end subroutine

!===========================================================================
! create image history
!  adj=F, input history file for reading data
!  adj=T, output history file for writing data
!---------------------------------------------------------------------------
  subroutine wei_imag_history(adj)
    logical, intent(in) :: adj
    integer             :: num_match
    character(len=1024) :: string
    !!.. create image history

    if (adj) then !! in case of migration 
       !!.. 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

       !!.. total number of extrapolation steps is set to be the same as the velocity model
       igrid%az   = vgrid%az

       call auxinout(img_)
       call auxputch("n1", "i", igrid%ax%n, img_)
       call auxputch("d1", "r", igrid%ax%d, img_)
       call auxputch("o1", "r", igrid%ax%o, img_)

       call auxputch("n2", "i", igrid%ay%n, img_)
       call auxputch("d2", "r", igrid%ay%d, img_)
       call auxputch("o2", "r", igrid%ay%o, img_)

       call auxputch("n3", "i", igrid%ahx%n, img_)
       call auxputch("d3", "r", igrid%ahx%d, img_)
       call auxputch("o3", "r", igrid%ahx%o, img_)

       call auxputch("n4", "i", igrid%ahy%n, img_)
       call auxputch("d4", "r", igrid%ahy%d, img_)
       call auxputch("o4", "r", igrid%ahy%o, img_)

       call auxputch("n5", "i", igrid%az%n, img_)
       call auxputch("d5", "r", igrid%az%d, img_)
       call auxputch("o5", "r", igrid%az%o, img_)

       call auxputch("esize" , "i", 4   , img_)
       call auxputch("label1", "s", "x" , img_)
       call auxputch("label2", "s", "y" , img_)
       call auxputch("label3", "s", "hx", img_)
       call auxputch("label4", "s", "hy", img_)
       call auxputch("label5", "s", "z" , img_)
    else ! in case of demigration
       call auxin(img_)
       call from_aux(img_, "n1", igrid%ax%n)
       call from_aux(img_, "d1", igrid%ax%d)
       call from_aux(img_, "o1", igrid%ax%o)
  
       call from_aux(img_, "n2", igrid%ay%n)
       call from_aux(img_, "d2", igrid%ay%d)
       call from_aux(img_, "o2", igrid%ay%o)

       call from_aux(img_, "n3", igrid%ahx%n)
       call from_aux(img_, "d3", igrid%ahx%d)
       call from_aux(img_, "o3", igrid%ahx%o)

       call from_aux(img_, "n4", igrid%ahy%n)
       call from_aux(img_, "d4", igrid%ahy%d)
       call from_aux(img_, "o4", igrid%ahy%o)

       call from_aux(img_, "n5", igrid%az%n)
       call from_aux(img_, "d5", igrid%az%d)
       call from_aux(img_, "o5", igrid%az%o)
    end if

    if (hessian_type > 0) then
       call auxinout(hessian_)
       call auxputch("n1", "i", igrid%ax%n, hessian_)
       call auxputch("d1", "r", igrid%ax%d, hessian_)
       call auxputch("o1", "r", igrid%ax%o, hessian_)

       call auxputch("n2", "i", igrid%ay%n, hessian_)
       call auxputch("d2", "r", igrid%ay%d, hessian_)
       call auxputch("o2", "r", igrid%ay%o, hessian_)

       call auxputch("n3", "i", igrid%ahx%n, hessian_)
       call auxputch("d3", "r", igrid%ahx%d, hessian_)
       call auxputch("o3", "r", igrid%ahx%o, hessian_)

       call auxputch("n4", "i", igrid%ahy%n, hessian_)
       call auxputch("d4", "r", igrid%ahy%d, hessian_)
       call auxputch("o4", "r", igrid%ahy%o, hessian_)

       call auxputch("n5", "i", igrid%az%n, hessian_)
       call auxputch("d5", "r", igrid%az%d, hessian_)
       call auxputch("o5", "r", igrid%az%o, hessian_)

       call auxputch("esize" , "i", 4   , hessian_)
       call auxputch("label1", "s", "x" , hessian_)
       call auxputch("label2", "s", "y" , hessian_)
       call auxputch("label3", "s", "hx", hessian_)
       call auxputch("label4", "s", "hy", hessian_)
       call auxputch("label5", "s", "z" , hessian_)

    end if

    !if (green_func) then
    !end if
  end subroutine 
!=====================================================================================
  subroutine wei_wfld_part_history(filetag, aws_in)
    character(len=*) :: filetag
    type(axis)         :: aws_in
    call auxputch("n1", "i", pgrid%ax%n, filetag)
    call auxputch("d1", "r", pgrid%ax%d, filetag)
    call auxputch("o1", "r", pgrid%ax%o, filetag)
 
    call auxputch("n2", "i", pgrid%ay%n, filetag)
    call auxputch("d2", "r", pgrid%ay%d, filetag)
    call auxputch("o2", "r", pgrid%ay%o, filetag)

    call auxputch("n3", "i", aws_in%n  , filetag)
    call auxputch("d3", "r", aws_in%d  , filetag)
    call auxputch("o3", "r", aws_in%o  , filetag)

    call auxputch("n4", "i", igrid%az%n, filetag)
    call auxputch("d4", "r", igrid%az%d, filetag)
    call auxputch("o4", "r", igrid%az%o, filetag)

    call auxputch("esize" , "i", 8     , filetag)
    call auxputch("label1", "s", "x"   , filetag)
    call auxputch("label2", "s", "y"   , filetag)
    call auxputch("label3", "s", "freq", filetag)
    call auxputch("label4", "s", "z"   , filetag)

  end subroutine
!=====================================================================================
  subroutine wei_wfld_full_history(filetag)
    character(len=*) :: filetag
    call auxputch("n1", "i", pgrid%ax%n, filetag)
    call auxputch("d1", "r", pgrid%ax%d, filetag)
    call auxputch("o1", "r", pgrid%ax%o, filetag)

    call auxputch("n2", "i", pgrid%ay%n, filetag)
    call auxputch("d2", "r", pgrid%ay%d, filetag)
    call auxputch("o2", "r", pgrid%ay%o, filetag)

    call auxputch("n3", "i", dgrid%aw%n, filetag)
    call auxputch("d3", "r", dgrid%aw%d, filetag)
    call auxputch("o3", "r", dgrid%aw%o, filetag)

    call auxputch("n4", "i", igrid%az%n, filetag)
    call auxputch("d4", "r", igrid%az%d, filetag)
    call auxputch("o4", "r", igrid%az%o, filetag)

    call auxputch("esize" , "i", 8     , filetag)
    call auxputch("label1", "s", "x"   , filetag)
    call auxputch("label2", "s", "y"   , filetag)
    call auxputch("label3", "s", "freq", filetag)
    call auxputch("label4", "s", "z"   , filetag)
  end subroutine
!=====================================================================================
  subroutine wei_io_file_init(adj)
    logical, intent(in) :: adj
    if (adj) then !!.. migration
       call wei_io_imag_init()
    else          !!.. demigration
       call wei_io_data_init()
    end if
    !! .. compute the diagonal of hessian
    if (hessian_type > 0) call wei_io_illu_init()
    !! .. Green's function
    !if (green_func) call wei_io_gren_init()
  end subroutine
!=====================================================================================
! image i/o initialization
!-------------------------------------------------------------------------------------
   subroutine wei_io_imag_init()
     real, allocatable, dimension(:,:,:,:) :: tmp
     integer                               :: iz, err
     allocate(tmp(igrid%ax%n, igrid%ay%n, igrid%ahx%n, igrid%ahy%n), stat=err)
     if (err /= 0) call seperr("trouble allocating tmp in wei_io_imag_init()")
     tmp = 0.
     do iz=1, vgrid%az%n
        call sep_write(tmp,img_, esize=4)
     end do
     deallocate(tmp)
   end subroutine
!=====================================================================================
   subroutine wei_io_illu_init()
     real, allocatable, dimension(:,:,:,:) :: tmp
     integer                               :: iz, err
     allocate(tmp(igrid%ax%n, igrid%ay%n, igrid%ahx%n, igrid%ahy%n), stat=err)
     if (err /= 0) call seperr("trouble allocating tmp in wei_io_imag_init()")
     tmp = 0.
     do iz=1, vgrid%az%n
        call sep_write(tmp,hessian_, esize=4)
     end do
     deallocate(tmp)
   end subroutine
!=====================================================================================
   subroutine wei_io_data_init()
     complex, allocatable, dimension(:,:,:) :: tmp
     integer                                :: ish, err
     allocate(tmp(dgrid%agx%n, dgrid%agy%n, dgrid%aw%n), stat=err)
     if (err /= 0) call seperr("trouble allocating tmp in wei_io_data_init()")
     tmp = 0.
     do ish=1, dgrid%nsh
        call sep_write(tmp, rec_, esize=8)
     end do
     deallocate(tmp)
   end subroutine
!=====================================================================================
! read or write one depth block image
!-------------------------------------------------------------------------------------
   subroutine wei_io_imag(flag, indz_in, nzs_in)
     character(len=4)      :: flag
     integer               :: indz_in
     integer               :: nzs_in
     integer, dimension(5) :: n_cube, j_wind, n_wind, f_wind
     integer               :: n_dim_cube, err
     integer               :: indz, n1

     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_in/)
     f_wind = (/ogrid%index_outimg_xmin-1, ogrid%index_outimg_ymin-1, 0, 0, indz_in-1/)


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

   end subroutine
!============================================================================
   subroutine wei_io_illu(flag, indz_in, nzs_in)
     character(len=4)      :: flag
     integer               :: indz_in
     integer               :: nzs_in
     integer, dimension(5) :: n_cube, j_wind, n_wind, f_wind
     integer               :: n_dim_cube, err
     integer               :: indz, n1

     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_in/)
     f_wind = (/ogrid%index_outimg_xmin-1, ogrid%index_outimg_ymin-1, 0, 0, indz_in-1/)

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

   end subroutine
!=============================================================================
!  read or write nws_in frequency slices of source wavefield 
!  for the indz_in'th depth step
!============================================================================
   subroutine wei_io_wfld(flag, filetag, wfld, nws_in, indz_in)
     character(len=4)      :: flag
     character(len=*)      :: filetag
     integer               :: nws_in, indw_in, indz_in
     integer, dimension(4) :: n_cube, j_wind, n_wind, f_wind
     integer               :: n_dim_cube, err
     complex, dimension(:,:,:) :: wfld
     n_dim_cube = 4
     n_cube = (/pgrid%ax%n, pgrid%ay%n, nws_in, pgrid%az%n/)
     j_wind = (/1, 1, 1, 1/)
     n_wind = (/pgrid%ax%n, pgrid%ay%n, nws_in, 1/)
     f_wind = (/0, 0, 0, indz_in-1/)

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

!=============================================================================
! read one frequency block receiver wavefield
!-----------------------------------------------------------------------------
   subroutine wei_io_crec(flag, indw_in, nws_in, ish_in)
     character(len=4)      :: flag
     integer               :: indw_in, nws_in, ish_in
     integer, dimension(3) :: n_cube, j_wind, n_wind, f_wind
     integer               :: n_dim_cube, err  
     integer               :: ishx, ishy

     n_dim_cube = 3
     n_cube = (/dgrid%ntr, dgrid%aw%n, dgrid%nsh/)
     j_wind = (/1, 1, 1/)
     n_wind = (/dgrid%ntr, nws_in, 1/)
     f_wind = (/0, indw_in-1, ish_in-1/)

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

!=============================================================================
! read one frequency block source wavefield
!-----------------------------------------------------------------------------
   subroutine wei_io_csou(indw_in, nws_in)
     integer               :: indw_in, nws_in
     integer, dimension(1) :: n_cube, j_wind, n_wind, f_wind
     integer               :: n_dim_cube, err
     n_dim_cube = 1
     n_cube = dgrid%aw%n
     j_wind = 1
     n_wind = nws_in
     f_wind = indw_in - 1
     csou = 0.
     err = sreed_window(sou_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,csou)
     if (err.lt.0) call seperr("wei_io_csou: trouble reading source wavefield")
   end subroutine

!==============================================================================
! read one depth block velocity model
!------------------------------------------------------------------------------
  subroutine wei_io_velo(indz_in, nzs_in)
    integer               :: indz_in
    integer               :: nzs_in
    integer, dimension(3) :: n_cube, j_wind, n_wind, f_wind
    integer               :: n_dim_cube, err, indz
    real, allocatable, dimension(:,:,:) :: tmp_vmod
    
    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_in/)
    f_wind = (/vgrid%index_provel_xmin-1, vgrid%index_provel_ymin-1, indz_in-1/)

    vmod = 0.
    err = sreed_window(vel_,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")

  end subroutine

!==============================================================================
  subroutine wei_io_head_init(flag)
    character(len=4), intent(in) :: flag
    if(flag == "hrec") then
       allocate(dgrid%header_rec(2, dgrid%ntr))
       allocate(pgrid%index_datpro_hrec(2, dgrid%ntr))
       dgrid%header_rec(2, dgrid%ntr) = 0.
       pgrid%index_datpro_hrec(2, dgrid%ntr) = 0.
    else if(flag == "hsou") then
       allocate(dgrid%header_sou(2))
       allocate(pgrid%index_datpro_hsou(2))
       dgrid%header_sou(2) = 0.
       pgrid%index_datpro_hsou(2) = 0.
    else
       call seperr("wei_io_head_init: you have to specify the type of header files")
    end if
  end subroutine
!==============================================================================
  subroutine wei_io_head_oper(flag)
    character(len=4), intent(in) :: flag
    if(flag == "hrec") then  
       call sep_read(dgrid%header_rec, hdr_rec_, esize=4)
    else if (flag == "hsou") then
       call sep_read(dgrid%header_sou, hdr_sou_, esize=4)
write(0,*)"in io_head_oper:", dgrid%header_sou(1), dgrid%header_sou(2)
    else
       call seperr("wei_io_head_lop: you have to specify the type of header files")
    end if
  end subroutine
!==============================================================================
!  subroutine wei_cwpsu_data(file_uid, index_trace)
!    integer             :: file_uid
!    integer             :: index_trace
!
!    real, dimension(60) :: header
!    integer(kind=4)     :: tracl             ! trace sequence number within line
!    integer(kind=4)     :: tracr             ! trace sequence number within reel
!    integer(kind=4)     :: fldr              ! field record number
!    integer(kind=4)     :: tracf             ! trace number within field record
!    integer(kind=4)     :: sx, sy, gx, gy
!    integer(kind=2)     :: ndt(2)            ! sample interval; in micro-seconds
!    integer(kind=2)     :: nlt(2)            ! number of samples in this trace

!    integer             :: itr, fldr_comp

!    equivalence(header(1) , tracl)
!    equivalence(header(2) , tracr)
!    equivalence(header(3) , fldr)
!    equivalence(header(4) , tracf)
!    equivalence(header(19), sx)
!    equivalence(header(20), sy)
!    equivalence(header(21), gx)
!    equivalence(header(22), gy)
!    equivalence(header(29), nlt(1))
!    equivalence(header(30), ndt(1))


!    header = 0.
!    itr = 0

    !!.. read the first header
!    read(file_uid, ref=index_trace)(header, ih=1, 60)
!    dgrid%sx = sx
!    dgrid%sy = sy
!    dgrid%nt = nlt(2)
!    dgrid%dt = ndt(1)/1000./1000.
!    fldr_comp = fldr
!    do while (fldr == fldr_comp)
!       read(file_uid, rec=index_trace+itr)(header, ih=1, 60)
!       itr = itr + 1
!    end do
!    ntr = itr

    !! have to compute the frequency: nw

!    do itr=1, ntr
!
!       read(file_uid, rec=index_trace+itr-1)(header, ih=1 , 60   )&
!                                            (ctrace, iws=1, aws%n)
!       crec_unsort(itr, :) = ctrace
!       header_rec(1, itr) = gx
!       header_rec(2, itr) = gy
!    end do

    

!  end subroutine
!==============================================================================
  

!=============================================================
! initiate file tags for input and output
!-------------------------------------------------------------
!   subroutine wei_io_tags_init()
!     tag_rec = "rec"
!     tag_sou = "sou"
!     tag_vel = "vel"
!     tag_img = "img"
!   end subroutine

end module










