module W_mod
  use down_mod
  use D_mod
  use R_mod
  !---------------
  implicit none

contains
  !----------------------------------------------------------------
  logical function bwfld_init(X,tag)
    type(twfld),       intent(inout) :: X
    character(len=128), intent(in)    :: tag

    bwfld_init=F
    X%tag = tag

    if(.not. wavefield_param_init(X,tag))then
       write(run%io,*) "trouble initializing wavefield by parameters"
       return
    end if
    
    bwfld_init=T
  end function bwfld_init
  !----------------------------------------------------------------
  logical function wavefield_history_init(X)
    type(twfld), intent(inout) :: X
#ifdef DEBUG
    character(len=128), parameter :: unit='wavefield_history_init'
    if(run%debug) call in(unit)
#endif
    wavefield_history_init=F
    call auxinout(X%tag)
    X%from = "history"

    call hinaxis(X%tag,X%amx,1); X%amx%l='mx'
    call hinaxis(X%tag,X%amy,2); X%amy%l='my'
    call hinaxis(X%tag,X%aw ,3); X%aw%l='w_'
    call hinaxis(X%tag,X%az ,4); X%az%l ='z_'
    call hinaxis(X%tag,X%asx,5); X%asx%l='sx'
    call hinaxis(X%tag,X%asy,6); X%asy%l ='sy'
    if(.not. wavefield_setup(X,.true.)) call seperr('wavefield_setup = notgood')
    wavefield_history_init=T
#ifdef DEBUG
    if(run%debug) call wavefield_report(X)
    if(run%debug) call out(unit)
#endif
  end function wavefield_history_init
  !----------------------------------------------------------------
  logical function create_wavefield_file(X,tag)
    type(twfld), intent(in) :: X
    integer:: er
    character(len=*)         :: tag
       create_wavefield_file=F 
       call auxscr(tag)
       call to_history("esize",8,tag)
       call houaxes(tag,X%amx,X%amy,X%aw,X%az,X%asx,X%asy)
       call param_write(tag)
       er=srite(X%tag,(0,0),8)
       er=sseek(X%tag,0,0)

       create_wavefield_file=T 
     end function create_wavefield_file
!----------------------------------------------------------------
   logical function wavefield_param_init(X,tag)
    type(twfld), intent(inout) :: X
    character(len=128), intent(in)    :: tag
#ifdef DEBUG
    character(len=128), parameter :: unit='wavefield_param_init'
    if(run%debug) call in(unit)
#endif
    wavefield_param_init=F

    X%from = "parameter"
    X%amx = down%amx;    X%amy = down%amy;
    X%asx = down%asx;    X%asy = down%asy;
    X%aw  = down%aw;
    X%az%d=down%az%d  
    X%az%n=down%az%n-down%firstzW+2 !+2 because includes z=0 and z=last
    X%az%o=(down%az%o) + (down%firstzW-1)*down%az%d 
    X%az%l='z_' ; X%aw%l='w_'; X%asy%l ='sy'
    X%amx%l='mx'; X%amy%l='my'; X%asx%l='sx'
    if(down%headD_d) X%asx=down%evnt
    if(.not. create_wavefield_file(X,X%tag)) then
       write(run%io,*) "trouble creating wavefield file"
       return
    end if
    if(.not.wavefield_setup(X,.true.)) then
      write(run%io,*) "trouble at wavefield_setup",X%tag
      return
    end if

    wavefield_param_init=T
    if(run%debug)call wavefield_report(X)
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end function wavefield_param_init
  !----------------------------------------------------------------
  logical function wavefield_setup(X,zero)
    type(twfld), intent(inout) :: X
    logical :: zero
    integer :: er
#ifdef DEBUG
    character(len=128), parameter :: unit='wavefield_setup'
    if(run%debug) call in(unit)
#endif
    wavefield_setup=F

    X%bites = 8* product( (/X%amx%n,X%amy%n,down%nws/) )

    call wavefield_report(X)
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
    wavefield_setup=T
  end function wavefield_setup

!----------------------------------------------------------------
  subroutine wavefield_assign(X,v)
    type(twfld)        :: X
    complex, intent(in) :: v
    logical ::a
#ifdef DEBUG
    character(len=128), parameter :: unit='wavefield_assign'
    if(run%debug) call in(unit)
#endif
    a=file_val_assign(X%tag,v)

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

    write(run%io,*)'        ************************************************'
    write(run%io,*)'        *       ',trim(X%tag),' from ',trim(X%from)
    write(run%io,*)'        *              [n]      [o]      [d]'
    write(run%io,'(a,i3,f10.2,f10.3)       ') "         *        Offx: ",X%amx%n, X%amx%o, X%amx%d
    write(run%io,'(a,i3,f10.2,f10.3)       ') "         *        Offy: ",X%amy%n, X%amy%o, X%amy%d
    write(run%io,'(a,i3,f10.2,f10.3)       ') "         *        Shtx: ",X%asx%n, X%asx%o, X%asx%d
    write(run%io,'(a,i3,f10.2,f10.3)       ') "         *        Shty: ",X%asy%n, X%asy%o, X%asy%d
    write(run%io,'(a,i3)       ') "         *        Total Shots : ",down%evnt%n
    write(run%io,'(a,i3,f10.2,f10.3,a,f5.2,a)')&
     "         *          W: ",X%aw%n,  X%aw%o,  X%aw%d," ->",real(down%aw%n)/real(down%nws)," blocks"
    write(run%io,'(a,i3,f10.2,f10.3,a,f5.2,a)')&
     "         *          Z: ",X%az%n,  X%az%o,  X%az%d," ->",real(down%az%n)/real(down%nzs)," blocks"

  end subroutine wavefield_report
!--------------------------------------------------------------
  logical function wavefield_read(X,data,iwb,izb,izs)
    type(twfld), intent(inout) :: X
    complex,dimension(:,:,:), intent(in)   :: data
    integer,intent(in) :: iwb,izb,izs
    integer :: stat,i1,i2,i3,ierr,n_dim_cube,fz
    integer,dimension(5)::n_cube,n_wind,f_wind,j_wind
#ifdef DEBUG
    character(len=128), parameter :: unit='wavefield_read'
    if(run%depar) call in(unit)
#endif
    wavefield_read=F
    fz=izs-1
    n_dim_cube=size(n_cube)
    n_cube=(/X%amx%n,X%amy%n,X%aw%n  ,X%az%n ,down%evnt%n/)
    j_wind=(/1      ,1      ,1       ,1      ,1           /)
    n_wind=(/X%amx%n,X%amy%n,down%nws,1      ,1           /)
    f_wind=(/0      ,0      ,(iwb-1)*down%wblock,(izb-1)*down%zblock+fz,down%shotnum /)
    if(run%lock) call sep_open_lock(X%tag)
    ierr=srite_window(X%tag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,data(1,1,1))
    if(run%lock) call sep_close_lock(X%tag)
    wavefield_read=T       
#ifdef DEBUG
    if(run%depar) call out(unit)
#endif
  end function wavefield_read
!----------------------------------------------------------------
!! writes out wavefield at before every extrapolation step
  logical function wavefield_write(X,data,iwb,izb,izs)
    type(twfld)                  :: X
    complex,dimension(:,:,:)    :: data
    integer,intent(in):: iwb,izb,izs
    integer :: stat,i1,i2,i3,ierr,n_dim_cube,fz
    integer,dimension(5)::n_cube,n_wind,f_wind,j_wind
#ifdef DEBUG
    character(len=128), parameter :: unit='wavefield_write'
    if(run%depar) call in(unit)
#endif
    wavefield_write=F
    fz=izs-1
    n_dim_cube=size(n_cube)
    n_cube=(/X%amx%n,X%amy%n,X%aw%n  ,X%az%n ,down%evnt%n /)
    j_wind=(/1      ,1      ,1       ,1      ,1            /)
    n_wind=(/X%amx%n,X%amy%n,down%nws,1      ,1            /)
    f_wind=(/0      ,0      ,(iwb-1)*down%wblock,(izb-1)*down%zblock+fz,down%shotnum /)

    if(run%lock) call sep_open_lock(X%tag)
    if (trim(X%tag)=='W') then
       ierr=srite_window(X%tag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,data(1,1,1))
       if(run%lock) call sep_close_lock(X%tag)
    elseif (trim(X%tag)=='W_d') then
       do i3=1,down%nws
          do i2=1,down%amy%n
             do i1=1,down%amx%n
                data(i1,i2,i3)=conjg(data(i1,i2,i3))
             end do
          end do
       end do
       ierr=srite_window(X%tag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,data(1,1,1))
       if(run%lock) call sep_close_lock(X%tag)
       do i3=1,down%nws
          do i2=1,down%amy%n
             do i1=1,down%amx%n
                data(i1,i2,i3)=conjg(data(i1,i2,i3))
             end do
          end do
       end do

    end if

    wavefield_write=T
#ifdef DEBUG
    if(run%depar) call out(unit)
#endif
  end function wavefield_write
!----------------------------------------------------------------
end module W_mod

