module wei_datum
  !! WEI migration

  use wei_util
  use wei_process
  use wei_work
   
  use wei_obj_S
  use wei_obj_R
  use wei_obj_D
  use wei_obj_W

  use wei_fft
  use wei_kmap
  use wei_xmap

  use wei_kxmig
  !---------------
  implicit none

contains
!!$
!!$=head1 
!!$
!!$ weidatum_init - Initialize datun operator
!!$
!!$=head1 SYNOPSIS
!!$
!!$  call weidatum_init(SLin,KXMIGin,FKin,FXin,IGin,in6,in7,in8,in9)
!!$
!!$=head1 PARAMETERS
!!$
!!$=over 4
!!$
!!$=item SLin - Slowness operator init  L<weislo1_init>
!!$
!!$=item KXMIGin - KXMIG operator init  L<weikxmig1_init>
!!$
!!$=item FKin - FK operator init  L<weiwem_init>
!!$
!!$=item FXin - FX operator init  L<weissf_init>
!!$
!!$=item IGin - image gather creation operator init  L<weihcig_init>
!!$
!!$=item in6,..in9 -  integer,optional
!!$
!!$ Dummy operators for future expansion
!!$
!!$=back
!!$
!!$=head1 DESCRIPTION
!!$
!!$ Initialize datuming operator
!!$
!!$=head1 SEE ALSO
!!$
!!$ L<weidatum> L<wei>
!!$
!!$=head1 LIBRARY
!!$
!!$ sepweif90
!!$
!!$=cut
!!$
  !----------------------------------------------------------------
  subroutine weidatum_init(SLin,KXMIGin,FKin,FXin,IGin,in6,in7,in8,in9)
    integer,external,optional ::                     in6,in7,in8,in9
    integer :: st,er
    character(len=128) :: name
#ifdef DEBUG
    character(len=128), parameter :: unit='weidatum_init'
#endif
    interface !!-------------------------
       integer function    SLin() result(st)
       end function SLin
       integer function KXMIGin() result(st)
       end function KXMIGin
       integer function    FKin() result(st)
       end function FKin
       integer function    FXin() result(st)
       end function FXin
       integer function    IGin() result(st)
       end function IGin
    end interface !!---------------------

    !! initialize the main WEI variables
    call process_init()
    
#ifdef DEBUG
    if(pro%debug) call in(unit)
#endif


    !! initialize objects
    call weibslow_init(S,pro%VS,pro%VG)  !INITIALIZE SLOW STRUCT WITH SLOWNESS TAG
    call distribute_init(.false.)!FIGURE OUT WHERE WE ARE GOING TO SEND STUFF
    call weibdata_init(D,pro%D)  !INITIALIZE DATA STRUCT WITH DATA TAG
    call weibdata_init(U,pro%W)  !INITIALIZE WAVEFIELD 

    !! auxiliary objects
#ifdef SEP_MPI
    write(name,"(a,a)") "bwfld.WEI",pro%tag_num
#else
    name="bwfld.WEI"  
#endif
    call weiauxdata_init(bW,name)   !CREATE AN AUXILARY WAVEFIELD FILE
    

    pro%mod_ = pro%D   !MODEL IS REFLECTIVITY
    pro%dat_ = pro%W   !DATA IS FREQUENCY SLICES

    !! convert frequency to radians
    call process_radians()

    !! init FFT
    call weifft_init()  

    !! init K-map
    call weikmap_init()
    call weixmap_init()


    !! allocate slowness arrays for FK/FX ops
    call weislow_allocate()

    !! check if offset is consistent with ph
    if(pro%pCIG) call check_hp()

    !! init operators
    st=SLin()     !! FK slownesses init
    st=KXMIGin()  !! KXMIG init
    st=FKin()     !! FK op init
    st=FXin()     !! FX op init
    st=IGin()     !! FX op init

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine weidatum_init
  !----------------------------------------------------------------
  !----------------------------------------------------------------
  subroutine weidatum_outnull()
  !ZERO OUTPUT IF REQUESTED
#ifdef DEBUG
    character(len=128), parameter :: unit='weidatum_outnull'
    if(pro%debug) call in(unit)
#endif

    if(pro%adj) then !! R<-D
       call      data_assign(U,(0.,0.))
    else             !! R->D
       call      data_assign(D,(0.,0.))
    end if

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine weidatum_outnull
  !----------------------------------------------------------------
!!$
!!$=head1 
!!$
!!$ weidatum -  Migration/modeling operator
!!$
!!$=head1 SYNOPSIS
!!$
!!$  ierr=weidatum(adj,add,mod,dat,SLop,KXMIGop,FKop,FXop,IGop,in6,in7,in8,in9)
!!$
!!$=head1 PARAMETERS
!!$
!!$=over 4
!!$
!!$=item adj,add - logcal
!!$
!!$ Adjoint and add (adjoint=.true. - down)
!!$ (adjoint=.false. - modeling).  add=.true. add to current file
!!$
!!$=item mod,dat - char*
!!$
!!$ Tags to the model and data
!!$
!!$=item SLop - Slowness operator   L<weislo1>
!!$
!!$=item KXMIGop - KXMIG operator   L<weikxmig1>
!!$
!!$=item FKop - FK operator   L<weiwem>
!!$
!!$=item FXop - FX operator   L<weissf>
!!$
!!$=item IGop - image gather creation operator   L<weihcig>
!!$
!!$=item in6,..in9 -  integer,optional
!!$
!!$ Dummy operators for future expansion
!!$
!!$=back
!!$
!!$=head1 DESCRIPTION
!!$
!!$ Migration/modeling operator
!!$
!!$ Three possible schemes  
!!$
!!$ Outer frequency---
!!$
!!$ iw_blocks{
!!$
!!$   iz_block{
!!$
!!$     migration
!!$
!!$   }
!!$
!!$ }
!!$
!!$ Outer depth---
!!$
!!$ iz_block{
!!$
!!$   iw_blocks{
!!$
!!$     migration
!!$
!!$   }
!!$
!!$ }
!!$
!!$ Outer frequency, variable steps---
!!$
!!$ iw_blocks{
!!$
!!$   iz_block (function of freq/slowness){
!!$
!!$     migration
!!$
!!$   }
!!$
!!$ }
!!$
!!$ If distribute specified 
!!$
!!$ distribute=vsf - variable steps as a function of frequency and slowness
!!$
!!$ distribute=vs - variable steps as a function of slowness
!!$
!!$ distribute=reg -  regular sampling methodology determination:
!!$
!!$ If outer_freq set then used otherwise
!!$
!!$ If few frequency blocks  outer frequency used
!!$
!!$ If few depth blocks  outer depth used
!!$
!!$
!!$=head1 SEE ALSO
!!$
!!$ L<weidatum_init>
!!$
!!$=head1 LIBRARY
!!$
!!$ sepweif90
!!$
!!$=cut
!!$
  !----------------------------------------------------------------
  integer function weidatum(adj,add,mod_,dat_,SLop,  KXMIGop,FKop,FXop,IGop,op6,op7,op8,op9)
  logical,            intent(in) :: adj,add
  character(len=128), intent(in) :: mod_,dat_
    integer,external,optional :: KXMIGop,FKop,FXop,IGop,op6,op7,op8,op9
    integer :: izb,iwb
    integer :: st
#ifdef DEBUG
    character(len=128), parameter :: unit='weidatum'
#endif
    interface !!----------------------------------------------
       integer function SLop() result(st)
       end function SLop
    end interface !!------------------------------------------
#ifdef DEBUG
    if(pro%debug) call in(unit)
#endif

    !! reset tag variables
    call process_reset(adj,add,mod_,dat_)
    if(pro%verb) call process_follow() !! where am I?
    
    !! accomodate tag name changes
#ifdef SEP_MPI
    R%tag_0 = pro%mod_
    D%tag_0 = pro%dat_
#else
    R%tag = trim(pro%mod_)
    D%tag = trim(pro%dat_)
#endif
   
#ifdef SEP_MPI
    if(pro%verb) write(0,*) "BEFORE DISTRIBUTE DATA",pro%mpi_ithread
    call distribute_data(D)
    if(pro%verb) write(0,*) "BEFORE DISTRIBUTE DOWN WAVE FIELD",pro%mpi_ithread
    call distribute_data(U)
    if(pro%verb) write(0,*) "AFTER DISTRIBUTE DOWN WAVE FIELD",pro%mpi_ithread
#endif


    !! zero the output if add is off
    if(.not. pro%add) call weidatum_outnull()

    !! setup the adjoint sign
    if(pro%adj) then
       pro%opsign = goT2B
    else
       pro%opsign = goB2T
    end if

    if(pro%outer_freq) then
       weidatum=weidatumof(SLop,  KXMIGop,FKop,FXop,IGop)
    else 
       weidatum=weidatumoz(SLop,  KXMIGop,FKop,FXop,IGop)
    end if

#ifdef SEP_MPI
    if(pro%verb) call weimessage("COMBINE DATA")
    call combine_data(D)
    if(pro%verb) call weimessage("OK")

    if(pro%verb) call weimessage("COMBINE DOWNWARD CONTINUATION")
    call combine_data(U)
    if(pro%verb) call weimessage("OK")
#endif

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
    weidatum=OK
    
  end function weidatum

  !----------------------------------------------------------
  integer function weidatumof(SLop,KXMIGop,FKop,FXop,IGop,op6,op7,op8,op9)
    integer,external,optional :: KXMIGop,FKop,FXop,IGop,op6,op7,op8,op9
    integer :: izb,iwb
    integer :: st
#ifdef DEBUG
    character(len=128), parameter :: unit='weidatumof'
#endif
    interface !!----------------------------------------------
       integer function SLop() result(st)
       end function SLop
    end interface !!------------------------------------------
     pro%dzb=pro%az%d

    if(pro%adj) then !! MIGRATION
#ifdef DEBUG
       if(pro%debug) call in('datuming')
#endif
       !----------------------------------------------------------
       call data_copy(U,bW)           !! copy the data at the top
       wbA: do iwb=1,pro%nwb
          call data_read(bW,iwb)      !! read  data at block top
          zbA: do izb=1,pro%nzb       
            call setup_var_block(iwb,izb)
            call slowness_read(S,izb) !! read a slowness block
            st = SLop()               !! setup sS, sR, s(z)
            if(pro%verb) call wei_followup(izb,iwb)
            R%bin=0.;
            pro%bw%o = pro%aw%o + (pro%iw_first(iwb)-1)*pro%bw%d
            call wdatum_down(izb,iwb,pro%nzs,KXMIGop,FKop,FXop,IGop)
            call data_write(bW,iwb) !! write data at block top
          end do zbA
       end do wbA
       call data_add(bW,D)          !! add U.C. wfld to data
       !----------------------------------------------------------
#ifdef DEBUG
       if(pro%debug) call out('datuming')
#endif
    else         !! MODELING
#ifdef DEBUG
       if(pro%debug) call in('modeling')
#endif
       !----------------------------------------------------------
       call data_copy(D,bW)           !! copy the data at the top
       wbF: do iwb=1,pro%nwb
           call data_read(bW,iwb)      !! read  data at block bottom
           zbF: do izb=pro%nzb,1,-1    
            call setup_var_block(iwb,izb)
             call    image_read(R,izb) !! read image block
             call slowness_read(S,izb) !! read a slowness block
             st = SLop()               !! setup sS,sR, s(z)
             if(pro%verb) call wei_followup(izb,iwb)
             pro%bw%o = pro%aw%o + (pro%iw_first(iwb)-1)*pro%bw%d
             call wdatum_down(izb,iwb,pro%nzs,KXMIGop,FKop,FXop,IGop) 
            end do zbF
            call data_write(bW,iwb) !! write data at block top
        end do wbF
       call data_add(bW,U)          !! add U.C. wfld to data
       !----------------------------------------------------------
#ifdef DEBUG
       if(pro%debug) call out('modeling')
#endif
    end if
#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
    weidatumof=OK
  end function weidatumof


  !----------------------------------------------------------
  integer function weidatumoz(SLop,KXMIGop,FKop,FXop,IGop,op6,op7,op8,op9)
    integer,external,optional :: KXMIGop,FKop,FXop,IGop,op6,op7,op8,op9
    integer :: izb,iwb
    integer :: st
#ifdef DEBUG
    character(len=128), parameter :: unit='weidatumoz'
#endif
    interface !!----------------------------------------------
       integer function SLop() result(st)
       end function SLop
    end interface !!------------------------------------------

       pro%dzb=pro%az%d
    if(pro%adj) then !! MIGRATION
#ifdef DEBUG
       if(pro%debug) call in('datum')
#endif
       !----------------------------------------------------------
       call data_copy(U,bW)         !! copy the data at the top
       zbA: do izb=1,pro%nzb        
          call slowness_read(S,izb) !! read a slowness block
          st = SLop()               !! setup sS, sR, s(z)
          wbA: do iwb=1,pro%nwb
            call setup_var_block(iwb,izb)
             if(pro%verb) call wei_followup(izb,iwb)
             pro%bw%o = pro%aw%o + (pro%iw_first(iwb)-1)*pro%bw%d
             call data_read(bW,iwb) !! read  data at block top
             call wdatum_down(izb,iwb,pro%nzs,KXMIGop,FKop,FXop,IGop) 
             call data_write(bW,iwb)  !! write data at block bottom
          end do wbA
       end do zbA
       call data_add(bw,D)     !! write image block
       !----------------------------------------------------------
#ifdef DEBUG
       if(pro%debug) call out('datuming')
#endif
    else         !! MODELING
#ifdef DEBUG
       if(pro%debug) call in('modeling')
#endif
       !----------------------------------------------------------
       call data_copy(D,bW)         !! copy the data at the top
       zbF: do izb=pro%nzb,1,-1
          call slowness_read(S,izb) !! read a slowness block
          st = SLop()               !! setup sS,sR, s(z)
          wbF: do iwb=1,pro%nwb
            call setup_var_block(iwb,izb)
             if(pro%verb) call wei_followup(izb,iwb)
             pro%bw%o = pro%aw%o + (pro%iw_first(iwb)-1)*pro%bw%d
             call data_read(bW,iwb) !! read  data at block bottom
                call wdatum_down(izb,iwb,pro%nzs,KXMIGop,FKop,FXop,IGop) 
             call data_write(bW,iwb)!! write data at block top
          end do wbF
       end do zbF
       call data_add(bW,U)          !! add U.C. wfld to data
       !----------------------------------------------------------
#ifdef DEBUG
       if(pro%debug) call out('modeling')
#endif
    end if
#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
    weidatumoz=OK
  end function weidatumoz
  !---------------------------------------------------------------
  subroutine wdatum_down(izb,iwb,nzs,KXMIGop,FKop,FXop,IGop)
    integer, intent(in)       :: izb,iwb
    integer,external,optional :: KXMIGop,FKop,FXop,IGop
    integer :: izs,iws,ith,nzs
#ifdef DEBUG
    character(len=128), parameter :: unit='wdatum_down'
    if(pro%debug) call in(unit)
#endif


    if(pro%adj) then !! MIGRATION
#ifdef SEP_OMP
       !$OMP  PARALLEL
       !$OMP  DO &
       !$OMP& SCHEDULE(DYNAMIC) &
       !$OMP& PRIVATE(iws,izs,ith)
#endif
       wsA: do iws=1,pro%nws             !LOOP OVER w STEPS
#ifdef SEP_OMP
          ith = omp_get_thread_num()+1
#else
          ith=1
#endif

          zsA: do izs=1,nzs          !LOOP OVER z STEPS

!             call weiimgW2R(izs,iws,ith,R%bin,bW%bin,IGop)
             call weikxmig(bW%bin,iws,izs,ith,KXMIGop,FKop,FXop)

          end do zsA
       end do wsA
#ifdef SEP_OMP
       !$OMP END DO
       !$OMP END PARALLEL
#endif
    else         !! MODELLING
#ifdef SEP_OMP
       !$OMP  PARALLEL
       !$OMP  DO &
       !$OMP& SCHEDULE(DYNAMIC) &
       !$OMP& PRIVATE(iws,izs,ith)
#endif
       wsF: do iws=1,pro%nws             !LOOP OVER w STEPS
#ifdef SEP_OMP
          ith = omp_get_thread_num()+1
#else
ith=1
#endif
          zsF: do izs=pro%nzs,1,-1       !LOOP OVER z STEPS
             call weikxmig(bW%bin,iws,izs,ith,KXMIGop,FKop,FXop)
!             call weiimgR2W(izs,iws,ith, R%bin,bW%bin,IGop)
          end do zsF
       end do wsF
#ifdef SEP_OMP
       !$OMP END DO
       !$OMP END PARALLEL
#endif
    end if

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


!!   ADJ         FOR
!!
!! D1 > R1     D1 > R1
!! v           ^
!! D2 > R2     D2 > R2
!! v           ^
!! D3 > R3     D3 > R3
!! v           ^
!! D4          D4
!! v           ^
!! ========    ========
!! v           ^
!! D4 > R4     D4 > R4
!! v           ^
!! D5 > R5     D5 > R5
!! v           ^
!! D6 > R6     D6 > R6
!! v           ^
!! D           D
!!



