module wei_obj_W
  use wei_util
  use wei_process
  use wei_type
  use wei_obj_D
  !---------------
  implicit none

  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!
  type(wewfld), public :: W
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!

contains
  !----------------------------------------------------------------
  subroutine weibwfld_init(X,tag)
    type(wewfld),       intent(inout) :: X
    character(len=128), intent(in)    :: tag

    select case(pro%operator(1:5))
    case('.Rw.D')

       X%tag = tag
       call wavefield_param_init(W)

    case('dSwdR')

       X%tag = tag
       if(isfile(X%tag)) then
          call wavefield_history_init(X)
       else
          call seperr('W /= file')
       end if
    end select

  end subroutine weibwfld_init
  !----------------------------------------------------------------
  subroutine wavefield_history_init(X)
    type(wewfld), intent(inout) :: X
#ifdef DEBUG
    character(len=128), parameter :: unit='wavefield_history_init'
    if(pro%debug) call in(unit)
#endif

    call auxinout(X%tag)
    X%from = "history"

    call hinaxis(X%tag,X%amx,1); X%amx%l='mx'
    if(pro%amx%n ==0) then
       pro%amx = X%amx
    else
       if(.not. coax( X%amx, pro%amx)) then
          call wavefield_report(X)
          call seperr('W.amx /= pro.amx')
       end if
    end if

    call hinaxis(X%tag,X%amy,2); X%amy%l='my'
    if(pro%amy%n ==0) then
       pro%amy = X%amy
    else
       if(.not. coax( X%amy, pro%amy)) then
          call wavefield_report(X)
          call seperr('W.amy /= pro.amy')
       end if
    end if

    call hinaxis(X%tag,X%ahx,3); X%ahx%l='hx'
    if(pro%ahx%n ==0) then
       pro%ahx = X%ahx
    else
       if(.not. coax( X%ahx, pro%ahx)) then
          call wavefield_report(X)
          call seperr('W.ahx /= pro.ahx')
       end if
    end if

    call hinaxis(X%tag,X%ahy,4); X%ahy%l='hy'
    if(pro%ahy%n ==0) then
       pro%ahy = X%ahy
    else
       if(.not. coax( X%ahy, pro%ahy)) then
          call wavefield_report(X)
          call seperr('W.ahy /= pro.ahy')
       end if
    end if

    call hinaxis(X%tag,X%aw ,5); X%aw%l ='w_'
    if(pro%aw%n ==0) then
       pro%aw = X%aw
    else
       if(.not. coax( X%aw, pro%aw)) then
          call wavefield_report(X)
          call seperr('W.aw /= pro.aw')
       end if
    end if

    call hinaxis(X%tag,X%az ,6); X%az%l ='z_'
    if(pro%az%n ==0) then
       pro%az = X%az
    else
       if(.not. coax( X%az, pro%az)) then
          call wavefield_report(X)
          call seperr('W.az /= pro.az')
       end if
    end if

    call wavefield_setup(X)
#ifdef DEBUG
    if(pro%debug) call wavefield_report(X)
    if(pro%debug) call out(unit)
#endif
  end subroutine wavefield_history_init
  !----------------------------------------------------------------
  subroutine wavefield_param_init(X)
    type(wewfld), intent(inout) :: X
#ifdef DEBUG
    character(len=128), parameter :: unit='wavefield_param_init'
    if(pro%debug) call in(unit)
#endif

    call auxscr(X%tag)
    X%from = "parameter"

    X%amx = pro%amx; call houaxis(X%tag,X%amx,1)
    X%amy = pro%amy; call houaxis(X%tag,X%amy,2)
    X%ahx = pro%ahx; call houaxis(X%tag,X%ahx,3)
    X%ahy = pro%ahy; call houaxis(X%tag,X%ahy,4)
    X%aw  = pro%aw;  call houaxis(X%tag,X%aw, 5)
    X%az  = pro%az;  call houaxis(X%tag,X%az, 6)

    call to_history("esize",8,X%tag)

    call wavefield_setup(X)
    call wavefield_assign(X,cmplx(0.))
#ifdef DEBUG
    if(pro%debug) call wavefield_report(X)
    if(pro%debug) call out(unit)
#endif
  end subroutine wavefield_param_init
  !----------------------------------------------------------------
  subroutine wavefield_setup(X)
    type(wewfld), intent(inout) :: X
    integer :: er
#ifdef DEBUG
    character(len=128), parameter :: unit='wavefield_setup'
    if(pro%debug) call in(unit)
#endif

    !! DEPTH blocking
    if( (1.0*X%az%n)/pro%nzs /= X%az%n/pro%nzs ) then
       call wavefield_report(W)
       call seperr('W.az.n is not a multiple of W.nzs')
    end if

    !! FREQUENCY blocking
    if( (1.0*X%aw%n)/pro%nws /= X%aw%n/pro%nws ) then
       call wavefield_report(W)
       call seperr('W.aw.n is not a multiple of W.nws')
    end if

    pro%nwb = pro%aw%n / pro%nws
    pro%nzb = pro%az%n / pro%nzs
    X%bites = 8* product( (/pro%amx%n,pro%amy%n,pro%ahx%n,pro%ahy%n/) )

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine wavefield_setup
  !----------------------------------------------------------------
  subroutine wavefield_assign(X,v)
    type(wewfld), intent(inout) :: X
    complex, intent(in) :: v
    integer :: izb, izs
    integer :: iwb, iws
#ifdef DEBUG
    character(len=128), parameter :: unit='wavefield_assign'
    if(pro%debug) call in(unit)
#endif

    bW%bin=v
    do izb=1,pro%nzb
       do izs=1,pro%nzs
          do iwb=1,pro%nwb
             do iws=1,pro%nws
                call wavefield_write(X,bW%bin,izb,izs,iwb,iws)
             end do
          end do
       end do
    end do

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine wavefield_assign
  !----------------------------------------------------------------  
  subroutine wavefield_report(X)
    type(wewfld), intent(in) :: X

    write(0,*)'        ************************************************'
    write(0,*)'        *       ',trim(X%tag),' from ',trim(X%from)
    write(0,'(a,i3,f10.2,f10.3)       ') "         *    Mx: ",X%amx%n, X%amx%o, X%amx%d
    write(0,'(a,i3,f10.2,f10.3)       ') "         *    My: ",X%amy%n, X%amy%o, X%amy%d
    write(0,'(a,i3,f10.2,f10.3)       ') "         *    Hx: ",X%ahx%n, X%ahx%o, X%ahx%d
    write(0,'(a,i3,f10.2,f10.3)       ') "         *    Hy: ",X%ahy%n, X%ahy%o, X%ahy%d
    write(0,'(a,i3,f10.2,f10.3,a,i3,a)') "         *     W: ",X%aw%n,  X%aw%o,  X%aw%d," ->",pro%nwb," blocks"
    write(0,'(a,i3,f10.2,f10.3,a,i3,a)') "         *     Z: ",X%az%n,  X%az%o,  X%az%d," ->",pro%nzb," blocks"
    write(0,*)'        ************************************************'
  end subroutine wavefield_report
  !--------------------------------------------------------------
  subroutine wavefield_read(X,data,izb,izs,iwb,iws)
    type(wewfld), intent(in) :: X
    complex, dimension(:,:,:,:,:), pointer :: data
    integer, intent(in) :: izb,izs,iwb,iws
    integer :: k
    integer :: stat
#ifdef DEBUG
    character(len=128), parameter :: unit='wavefield_read'
    if(pro%depar) call in(unit)
#endif

    k=(izb-1)*pro%nzs*pro%nwb*pro%nws &
    +         (izs-1)*pro%nwb*pro%nws &
    +                 (iwb-1)*pro%nws &
    +                             iws

    stat=sseek_block(X%tag,k-1,X%bites,0)
    call sreed(X%tag, data(:,:,:,:,iws), X%bites)

#ifdef DEBUG
    if(pro%depar) call out(unit)
#endif
  end subroutine wavefield_read
  !----------------------------------------------------------------
  subroutine wavefield_write(X,data,izb,izs,iwb,iws)
    type(wewfld), intent(in) :: X
    complex, dimension(:,:,:,:,:), pointer :: data
    integer, intent(in) :: izb,izs,iwb,iws
    integer :: k
    integer :: stat
#ifdef DEBUG
    character(len=128), parameter :: unit='wavefield_write'
    if(pro%depar) call in(unit)
#endif

    k=(izb-1)*pro%nzs*pro%nwb*pro%nws &
    +         (izs-1)*pro%nwb*pro%nws &
    +                 (iwb-1)*pro%nws &
    +                             iws

    stat=sseek_block(X%tag,k-1,X%bites,0)
    call srite(X%tag, data(:,:,:,:,iws), X%bites)

#ifdef DEBUG
    if(pro%depar) call out(unit)
#endif
  end subroutine wavefield_write
  !----------------------------------------------------------------
end module wei_obj_W
