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

  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!
  type(wedata), public :: D      ! background   data
  type(wedata), public :: U      ! Datumed wavefield
  type(wedata), public ::bW      ! background   wavefield
  type(wedata), public ::pW      ! perturbation wavefield
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!

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

    X%tag = tag
    if(isfile(X%tag)) then
       call data_history_init(X)
    else
       call data_param_init(X)
    end if

  end subroutine weibdata_init
  !----------------------------------------------------------------
  subroutine weiauxdata_init(X,tag)
    type(wedata),       intent(inout) :: X
    character(len=128), intent(in)    :: tag

    X%tag = tag
    call data_scratch_init(X)

  end subroutine weiauxdata_init
  !----------------------------------------------------------------
  subroutine data_history_init(X)
    type(wedata), intent(inout) :: X
#ifdef DEBUG
    character(len=128), parameter :: unit='data_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 data_report(X)
          call seperr('D.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 data_report(X)
          call seperr('D.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 data_report(X)
          call seperr('D.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 data_report(X)
          call seperr('D.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 data_report(X)
          call seperr('D.aw /= pro.aw')
       end if
    end if

    X%aw_local=X%aw; X%aw_local%n=pro%nw_local
    call data_setup(X)

#ifdef SEP_MPI
    if(pro%mpi_nthread>1) then
       X%tag_0=X%tag;
       write(X%tag,"(a,a)") trim(X%tag_0),pro%tag_num
       call auxscr(X%tag)
       call houaxis(X%tag,X%amx,1);call houaxis(X%tag,X%amy,2)
       call houaxis(X%tag,X%ahx,3);call houaxis(X%tag,X%ahy,4)
       call houaxis(X%tag,X%aw_local, 5)
       call to_history("esize",8,X%tag);
    end if
#endif

#ifdef DEBUG
    if(pro%debug) call data_report(X)
    if(pro%debug) call out(unit)
#endif
  end subroutine data_history_init
  !----------------------------------------------------------------
  subroutine data_param_init(X)
    type(wedata), intent(inout) :: X
#ifdef DEBUG
    character(len=128), parameter :: unit='data_param_init' 
    if(pro%debug) call in(unit)
#endif

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

    X%amx = pro%amx; X%amy = pro%amy; 
    X%ahx = pro%ahx; X%ahy = pro%ahy; 
    X%aw  = pro%aw;  

    X%aw_local=X%aw; X%aw_local%n=pro%nw_local

#ifdef SEP_MPI
    if(pro%mpi_nthread>1) then
       X%tag_0=X%tag;
       write(X%tag,"(a,a)") trim(X%tag_0),pro%tag_num
       call auxscr(X%tag)
       call houaxis(X%tag_0,X%amx,1);call houaxis(X%tag_0,X%amy,2)
       call houaxis(X%tag_0,X%ahx,3);call houaxis(X%tag_0,X%ahy,4)
       call to_history("esize",8,X%tag_0)
       call houaxis(X%tag_0,X%aw, 5)
    end if
#endif

    call to_history("esize",8,X%tag)
    call houaxis(X%tag,X%amx,1);call houaxis(X%tag,X%amy,2)
    call houaxis(X%tag,X%ahx,3);call houaxis(X%tag,X%ahy,4)
    call houaxis(X%tag,X%aw_local,5)

    call data_setup(X)



    !THIS IS GOING TO BE A PROBLEM!!!
#ifdef SEP_MPI
    if(pro%mpi_ithread==0) call data_mpi_assign(X,cmplx(0.))
#else
    call data_assign(X,cmplx(0.))
#endif

#ifdef DEBUG
    if(pro%debug) call data_report(X)
    if(pro%debug) call out(unit)
#endif
  end subroutine data_param_init
  !----------------------------------------------------------------
  subroutine data_scratch_init(X)
    type(wedata), intent(inout) :: X
#ifdef DEBUG
    character(len=128), parameter :: unit='data_scratch_init' 
    if(pro%debug) call in(unit)
#endif

    call auxtmp(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;  
    X%aw_local=X%aw; X%aw_local%n=pro%nw_local
    call houaxis(X%tag,X%aw_local, 5)

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

    call data_setup(X)
    call data_assign(X,(0.,0.))
#ifdef DEBUG
    if(pro%debug) call data_report(X)
    if(pro%debug) call out(unit)
#endif
  end subroutine data_scratch_init
  !----------------------------------------------------------------
  subroutine data_setup(X)
    type(wedata), intent(inout) :: X
    integer :: er
#ifdef DEBUG
    character(len=128), parameter :: unit='data_setup' 
    if(pro%debug) call in(unit)
#endif

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

    pro%nwb = pro%nw_local / pro%nws
    !!--------------------
    allocate(X%bin(X%amx%n,pro%amy%n,pro%ahx%n,pro%ahy%n,pro%nws),stat=er)
    if(er/=0) call seperr('cant allocate D.bin')
    X%bin=(0.0,0.0)
    !!--------------------
    X%bites = 8* product( (/pro%amx%n,pro%amy%n,pro%ahx%n,pro%ahy%n,pro%nws/) )

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine data_setup
  !----------------------------------------------------------------
  subroutine data_mpi_assign(X,v)
    type(wedata), intent(inout) :: X
    complex, intent(in) :: v
    integer :: iwb,stat

    X%bin=v
    do iwb=1,size(pro%iw_thread)
       stat=sseek_block(X%tag_0,iwb-1,X%bites,0)
       call srite(X%tag_0, X%bin, X%bites)
    end do
   call sseek(X%tag_0,0,0)

  end subroutine data_mpi_assign
  !----------------------------------------------------------------
  subroutine data_assign(X,v)
    type(wedata), intent(inout) :: X
    complex, intent(in) :: v
    integer :: iwb
#ifdef DEBUG
    character(len=128), parameter :: unit='data_assign' 
    if(pro%debug) call in(unit)
#endif

    X%bin=v
    do iwb=1,pro%nwb
       call data_write(X,iwb)
    end do

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine data_assign
  !----------------------------------------------------------------
  subroutine data_report(X)
    type(wedata), intent(in) :: X
    write(0,*)'        ************************************************'
    write(0,*)'        *       ',trim(X%tag),' from ',trim(X%from)
    write(0,'(a,i3,f14.2,f10.3)       ') "         *    Mx: ",X%amx%n, X%amx%o, X%amx%d
    write(0,'(a,i3,f14.2,f10.3)       ') "         *    My: ",X%amy%n, X%amy%o, X%amy%d
    write(0,'(a,i3,f14.2,f10.3)       ') "         *    Hx: ",X%ahx%n, X%ahx%o, X%ahx%d
    write(0,'(a,i3,f14.2,f10.3)       ') "         *    Hy: ",X%ahy%n, X%ahy%o, X%ahy%d
    write(0,'(a,i3,a,i3,a,f10.2,f10.3,a,i3,a)') "         *     W: ",X%aw%n,&
    "(",X%aw_local%n,")",  X%aw%o,  X%aw%d," ->",pro%nwb," blocks"
    if(associated(X%bin)) then
       write(0,'(a,f9.2,a,f9.2,a)        ') "         * block: ",X%bites/1024.," / ",8*size(X%bin)/1024.," Kbytes"
    end if
    write(0,*)'        ************************************************'
  end subroutine data_report
  !--------------------------------------------------------------
  subroutine data_read(X,k)
    type(wedata), intent(inout) :: X
    integer, intent(in) :: k
    integer :: stat
#ifdef DEBUG
    character(len=128), parameter :: unit='data_read'
    if(pro%debug) call in(unit)
#endif

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

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine data_read
  !----------------------------------------------------------------
  subroutine data_write(X,k)
    type(wedata), intent(inout) :: X
    integer, intent(in) :: k
    integer :: stat
#ifdef DEBUG
    character(len=128), parameter :: unit='data_write' 
    if(pro%debug) call in(unit)
#endif

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

     call to_history("ifreq",k,X%tag)

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine data_write
  !----------------------------------------------------------------
  subroutine distribute_data(X)
    type(wedata), intent(inout) :: X
#ifdef SEP_MPI
    if(pro%mpi_nthread>1) then
       if(0/=mpi_sep_tag_distribute(0,X%tag_0,pro%mpi_ithread,X%tag,X%bites/4,&
       pro%nwp, pro%iw_thread)) &
       call seperr("trouble distributing data ");
    end if
#endif
  end subroutine distribute_data
  !----------------------------------------------------------------
  subroutine combine_data(X)
    type(wedata), intent(inout) :: X
    integer :: ierr,i
    integer, external :: sseek
#ifdef SEP_MPI
     i=0;if(pro%add) i=1
    if( .not. pro%adj .and. pro%mpi_nthread>1) then
       if(pro%mpi_ithread==0) ierr= sseek(trim(x%tag_0),0,0)
       ierr=sseek(trim(X%tag),0,0)

       if(0/=mpi_sep_tag_combine(pro%mpi_ithread,X%tag,0,X%tag_0, &
       X%bites/4,pro%nwp,pro%iw_thread,i)) &
       call seperr("trouble combining data ");
    end if
    call aux_unlink(X%tag)
#endif
  end subroutine combine_data
  !----------------------------------------------------------------
  subroutine data_copy(X,Y)
    type(wedata), intent(inout) :: X
    type(wedata), intent(inout) :: Y
    integer :: iwb
#ifdef DEBUG
    character(len=128), parameter :: unit='data_copy'
    if(pro%debug) call in(unit)
#endif

    do iwb=1,pro%nwb
       call data_read (X,iwb)
       Y%bin = X%bin
       call data_write(Y,iwb)
    end do

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine data_copy
  !----------------------------------------------------------------
  subroutine data_add(X,Y)
    type(wedata), intent(inout) :: X
    type(wedata), intent(inout) :: Y
    integer :: iwb
#ifdef DEBUG
    character(len=128), parameter :: unit='data_add'
    if(pro%debug) call in(unit)
#endif

    do iwb=1,pro%nwb
       call data_read (X,iwb)
       call data_read (Y,iwb)
       Y%bin = Y%bin + X%bin
       call data_write(Y,iwb)
    end do

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine data_add
  !----------------------------------------------------------------  
end module wei_obj_D









