module wei_mig
  !! 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
  use wei_image
  !---------------
  implicit none

  logical, private :: wIO
contains
!!$
!!$=head1 
!!$
!!$ weimig_init - Initialize migration operator
!!$
!!$=head1 SYNOPSIS
!!$
!!$  call weimig_init(SLin,WCin,FKin,FXin,IGin,in6,in7,in8,in9)
!!$
!!$=head1 PARAMETERS
!!$
!!$=over 4
!!$
!!$=item SLin - Slowness operator init  L<weislo1_init>
!!$
!!$=item WCin - WCop 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 migration/modeling operator
!!$
!!$=head1 SEE ALSO
!!$
!!$ L<weimig> L<wei>
!!$
!!$=head1 LIBRARY
!!$
!!$ sepweif90
!!$
!!$=cut
!!$
  !----------------------------------------------------------------
  subroutine weimig_init(SLin,WCin,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='weimig_init'
#endif
    interface !!-------------------------
       integer function    SLin() result(st)
       end function SLin
       integer function WCin() result(st)
       end function WCin
       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

    !! flag whether I/O on wavefield.
    wIO=.false.
    select case(pro%operator(1:5))
    case('.Rw.D')
       wIO=.true.
    end select

    !! initialize objects
write(0,*) "I DIE 1"
    call weibslow_init(S,pro%VS,pro%VG)!INITIALIZE SLOW STRUCT WITH SLOWNESS TAG
write(0,*) "I DIE 2"
    call distribute_init()       !FIGURE OUT WHERE WE ARE GOING TO SEND STUFF
write(0,*) "I DIE 3"
    call weibimag_init(R,pro%R)  !INITIALIZE REFL. STRUCT WITH REFLECTIVITY TAG
write(0,*) "I DIE 4"
    call weibdata_init(D,pro%D)  !INITIALIZE DATA STRUCT WITH DATA TAG
write(0,*) "I DIE 5"

    !! 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
    
    if(wIO) call weibwfld_init(W,pro%W) !! needs aux obj  

    pro%mod_ = pro%R   !MODEL IS REFLECTIVITY
    pro%dat_ = pro%D   !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()

    !! init imaging
    call weiimage_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=WCin()  !! WCop 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 weimig_init
  !----------------------------------------------------------------
  !----------------------------------------------------------------
  subroutine weimig_outnull()
  !ZERO OUTPUT IF REQUESTED
#ifdef DEBUG
    character(len=128), parameter :: unit='weimig_outnull'
    if(pro%debug) call in(unit)
#endif

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

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine weimig_outnull
  !----------------------------------------------------------------
!!$
!!$=head1 
!!$
!!$ weimig -  Migration/modeling operator
!!$
!!$=head1 SYNOPSIS
!!$
!!$  ierr=weimig(adj,add,mod,dat,SLop,WCop,FKop,FXop,IGop,in6,in7,in8,in9)
!!$
!!$=head1 PARAMETERS
!!$
!!$=over 4
!!$
!!$=item adj,add - logcal
!!$
!!$ Adjoint and add (adjoint=.true. - migration)
!!$ (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 WCop - WCop 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<weimig_init>
!!$
!!$=head1 LIBRARY
!!$
!!$ sepweif90
!!$
!!$=cut
!!$
  !----------------------------------------------------------------
  integer function weimig(adj,add,mod_,dat_,SLop,  WCop,FKop,FXop,IGop,op6,op7,op8,op9)
  logical,            intent(in) :: adj,add
  character(len=128), intent(in) :: mod_,dat_
    integer,external,optional :: WCop,FKop,FXop,IGop,op6,op7,op8,op9
    integer :: izb,iwb
    integer :: st
#ifdef DEBUG
    character(len=128), parameter :: unit='weimig'
#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(pro%io_unit,*) "BEFORE DISTRIBUTE DATA",pro%mpi_ithread
    call distribute_data(D)
    if(pro%verb) write(pro%io_unit,*) "BEFORE DISTRIBUTE IMAGE",pro%mpi_ithread
    call distribute_image(R)
    if(pro%verb) write(pro%io_unit,*) "AFTER DISTRIBUTE IMAGE",pro%mpi_ithread
#endif


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

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

    call to_history("nz",pro%nzb,R%tag) 
    call to_history("nz",pro%nzb,R%tag) 
    if(pro%distribute_pattern>0) then
       weimig=weimigvar(SLop,  WCop,FKop,FXop,IGop)
    else if(pro%outer_freq) then
       weimig=weimigof(SLop,  WCop,FKop,FXop,IGop)
    else 
       weimig=weimigoz(SLop,  WCop,FKop,FXop,IGop)
    end if

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

    if(pro%verb) call weimessage("COMBINE IMAGE")
    call combine_image(R)

    if(pro%verb) call weimessage("COMBINE SLOW")
    call combine_slow(S)

    if(pro%verb) call weimessage("OK")
#endif

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

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

    if(pro%adj) then !! MIGRATION
#ifdef DEBUG
       if(pro%debug) call in('migration')
#endif
       !----------------------------------------------------------
       call data_copy(D,bW)           !! copy the data at the top
          call to_history("fblock",0,R%tag) 
       wbA: do iwb=1,pro%nwb
          call data_read(bW,iwb)      !! read  data at block top
          call var_samp_iw(iwb)
          call var_slow(iwb)
          izdone=0;izb=0
          nzb=ceiling(real(size(pro%z_samp(iwb)%dz))/real(pro%nzs))
          zbA: do while(izdone< size(pro%z_samp(iwb)%dz)) 
            izb=izb+1
            call setup_var_block(iwb,izb)
            nzs=min(pro%nzs,size(pro%z_samp(iwb)%dz)-izdone)
            R%bin=cmplx(0.,0.)
            call slowness_read(S,izb,nzs) !! read a slowness block
            st = SLop()               !! setup sS, sR, s(z)
            if(pro%verb) call wei_followup(izb,iwb,nzb)

            pro%bw%o = pro%aw%o + (pro%iw_first(iwb)-1)*pro%bw%d
            call wmigRD(izb,iwb,nzs,WCop,FKop,FXop,IGop) !! trash W
            call image_write(R,izb,iwb/=1.or.pro%add,iwb,nzs)!!write image block
            izdone=izdone+nzs

          end do zbA
          call to_history("fblock",iwb,R%tag) 
       end do wbA
       !----------------------------------------------------------
#ifdef DEBUG
       if(pro%debug) call out('migration')
#endif
    else         !! MODELING

#ifdef DEBUG
       if(pro%debug) call in('modeling')
#endif
       call seperr("unimplemented variable sampling modeling")
#ifdef DEBUG
       if(pro%debug) call out('modeling')
#endif
    end if
#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
    weimigvar=OK
  end function weimigvar
  !----------------------------------------------------------
  integer function weimigof(SLop,WCop,FKop,FXop,IGop,op6,op7,op8,op9)
    integer,external,optional :: WCop,FKop,FXop,IGop,op6,op7,op8,op9
    integer :: izb,iwb
    integer :: st
#ifdef DEBUG
    character(len=128), parameter :: unit='weimigof'
#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('migration')
#endif
       !----------------------------------------------------------
       call data_copy(D,bW)           !! copy the data at the top
       call to_history("fblock",0,R%tag) 
       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
            if(wIO) then           !! DC one block
                call wmigWW(izb,iwb,WCop,FKop,FXop,IGop) !! write W
             else
                call wmigRD(izb,iwb,pro%nzs,WCop,FKop,FXop,IGop) !! trash W
             end if
            call image_write(R,izb,(iwb/=1 .or. pro%add)) !! write image block
          end do zbA
          call to_history("fblock",iwb,R%tag) 
       end do wbA
       !----------------------------------------------------------
#ifdef DEBUG
       if(pro%debug) call out('migration')
#endif
    else         !! MODELING
#ifdef DEBUG
       if(pro%debug) call in('modeling')
#endif
       !----------------------------------------------------------
       call data_assign(bW,(0.,0.))    !! start with no data
       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 + (iwb-1)*(pro%bw%n*pro%bw%d)
             pro%bw%o = pro%aw%o + (pro%iw_first(iwb)-1)*pro%bw%d
             if(wIO) then           !! DC one block
                call wmigWW(izb,iwb,WCop,FKop,FXop,IGop) !! write W
             else
                call wmigRD(izb,iwb,pro%nzs,WCop,FKop,FXop,IGop) !! trash W
             end if
            end do zbF
            call data_write(bW,iwb) !! write data at block top
             call to_history("fblock",iwb,D%tag) 
        end do wbF
       call data_add(bW,D)          !! 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
    weimigof=OK
  end function weimigof


  !----------------------------------------------------------
  integer function weimigoz(SLop,WCop,FKop,FXop,IGop,op6,op7,op8,op9)
    integer,external,optional :: WCop,FKop,FXop,IGop,op6,op7,op8,op9
    integer :: izb,iwb
    integer :: st
#ifdef DEBUG
    character(len=128), parameter :: unit='weimigoz'
#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('migration')
#endif
       !----------------------------------------------------------
       call data_copy(D,bW)         !! copy the data at the top
       zbA: do izb=1,pro%nzb        
          call    image_read(R,izb) !! read image block
          call slowness_read(S,izb) !! read a slowness block
          st = SLop()               !! setup sS, sR, s(z)
          call to_history("fblock",0,R%tag) 
          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
             if(wIO) then           !! DC one block
                call wmigWW(izb,iwb,WCop,FKop,FXop,IGop) !! write W
             else
                call wmigRD(izb,iwb,pro%nzs,WCop,FKop,FXop,IGop) !! trash W
             end if
             call data_write(bW,iwb)  !! write data at block bottom
             call to_history("fblock",iwb,R%tag) 
          end do wbA
          call image_write(R,izb)     !! write image block
       end do zbA
       !----------------------------------------------------------
#ifdef DEBUG
       if(pro%debug) call out('migration')
#endif
    else         !! MODELING
#ifdef DEBUG
       if(pro%debug) call in('modeling')
#endif
       !----------------------------------------------------------
       call data_assign(bW,(0.,0.)) !! start with no data
       zbF: do izb=pro%nzb,1,-1
          call    image_read(R,izb) !! read image block
          call slowness_read(S,izb) !! read a slowness block
          st = SLop()               !! setup sS,sR, s(z)
          call to_history("fblock",0,D%tag) 
          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 + (iwb-1)*(pro%bw%n*pro%bw%d)
             pro%bw%o = pro%aw%o + (pro%iw_first(iwb)-1)*pro%bw%d
             call data_read(bW,iwb) !! read  data at block bottom
             if(wIO) then           !! DC one block
                call wmigWW(izb,iwb,WCop,FKop,FXop,IGop) !! write W
             else
                call wmigRD(izb,iwb,pro%nzs,WCop,FKop,FXop,IGop) !! trash W
             end if
             call data_write(bW,iwb)!! write data at block top
             call to_history("fblock",iwb,D%tag) 
          end do wbF
       end do zbF
       call data_add(bW,D)          !! 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
    weimigoz=OK
  end function weimigoz
  !---------------------------------------------------------------
  !! don't write the wavefield
  subroutine wmigRD(izb,iwb,nzs,WCop,FKop,FXop,IGop)
    integer, intent(in)       :: izb,iwb
    integer,external,optional :: WCop,FKop,FXop,IGop
    integer :: izs,iws,ith,nzs
#ifdef DEBUG
    character(len=128), parameter :: unit='wmigRD'
    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,WCop,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,WCop,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 wmigRD
  !----------------------------------------------------------------
  !! write the wavefield
  subroutine wmigWW(izb,iwb,WCop,FKop,FXop,IGop)
    integer, intent(in)       :: izb,iwb
    integer,external,optional :: WCop,FKop,FXop,IGop
    integer :: izs,iws,ith
#ifdef DEBUG
    character(len=128), parameter :: unit='wmigWW'
    if(pro%debug) call in(unit)
#endif

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

             !$OMP CRITICAL
             call wavefield_write(W,bW%bin,izb,izs,iwb,iws)
             !$OMP END CRITICAL

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

          end do zsA
       end do wsA
       !$OMP END DO
       !$OMP END PARALLEL
    else         !! MODELLING
       !$OMP  PARALLEL
       !$OMP  DO &
       !$OMP& SCHEDULE(DYNAMIC) &
       !$OMP& PRIVATE(iws,izs,ith)
       wsF: do iws=1,pro%nws
#ifdef SEP_OMP
          ith = omp_get_thread_num()+1
#else
          ith = 1
#endif
          zsF: do izs=pro%nzs,1,-1

             call weikxmig(bW%bin,iws,izs,ith,WCop,FKop,FXop)
             call weiimgR2W(izs,iws,ith, R%bin,bW%bin,IGop)

             !$OMP CRITICAL
             call wavefield_write(W,bW%bin,izb,izs,iwb,iws)
             !$OMP END CRITICAL

          end do zsF
       end do wsF
       !$OMP END DO
       !$OMP END PARALLEL
    end if

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


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



