module wei_dtm
  !! WEI datuming

  use wei_util
  use wei_pro
  use wei_par
  use wei_tag
  use wei_comp

  use wei_S
  use wei_D

  use wei_fft
  use wei_kmap
  use wei_xmap

  use wei_kxmig

  use wei_mig


  !---------------
  implicit none

  type(weslow), private :: SS
  type(wedata), private :: DD 
  type(wedata), private :: UU

  type(wedata), private :: cW



contains
  !----------------------------------------------------------------
  subroutine weidtm_init(SLin,WCin,FKin,FXin,in5,in6,in7,in8,in9)
    integer,external,optional ::             in5,in6,in7,in8,in9
    integer :: st,er
    character(len=128) :: name
#ifdef DEBUG
    character(len=128), parameter :: unit='weidtm_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
    end interface !!---------------------

#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    wpro%mod_ = wtag%D
    wpro%dat_ = wtag%U

    cW%tag="cW"; call weiD_aux_init(cW)

    !! initialize objects
    SS%tag=wtag%S; call weiS_init(SS)
    DD%tag=wtag%D; call weiD_init(DD)
    UU%tag=wtag%U; call weiD_init(UU)

    !! convert frequency to radians
    call weipro_radians()

    call weifft_init()    !! init FFT
    call weikmap_init()   !! init K-map
    call weixmap_init()   !! init X-map

    !! allocate slowness arrays for FK/FX ops
    call weiS_fkfx_allocate()
    
    !! init source/receiver pointers
    call weicomp_init()

    !! check if offset is consistent with ph
    if(wpro%pCIG) call weipro_check_hp()

    !! init operators
    st=SLin()     !! SL op init
    st=WCin()     !! WC op init
    st=FKin()     !! FK op init
    st=FXin()     !! FX op init

    !! spread slowness on the nodes
    call weiS_mpicopy(SS)

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weidtm_init
  !----------------------------------------------------------------
  !----------------------------------------------------------------
  subroutine weidtm_outnull()
#ifdef DEBUG
    character(len=128), parameter :: unit='weidtm_outnull'
    if(wpro%dbg>0) call in(unit)
#endif

    if(wpro%adj) then !! D<-U
       call weiD_assign(DD,(0.,0.))
    else              !! D->U
       call weiD_assign(UU,(0.,0.))
    end if

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weidtm_outnull
  !----------------------------------------------------------------
  !----------------------------------------------------------------
  integer function weidtm(adj,add,mod_,dat_,SLop,WCop,FKop,FXop,op5,op6,op7,op8,op9)
    integer,external,optional      ::       SLop,WCop,FKop,FXop,op5,op6,op7,op8,op9
    logical,            intent(in) :: adj,add
    character(len=128), intent(in) :: mod_,dat_
    integer :: izb,iwb
    integer :: st
#ifdef DEBUG
    character(len=128), parameter :: unit='weidtm'
    if(wpro%dbg>0) call in(unit)
#endif

    !! reset tag variables
    call weipro_reset(adj,add,mod_,dat_)
    if(wpar%impi==0 .and. wpro%verb) call weipro_follow() !! where am I?

    !! accomodate tag name changes
    !! this is only needed on the master
    !! since the files on the nodes are temporary
    !! and going to be filled next
    if(wpar%impi==0) then
       DD%tag = trim(wpro%mod_)
       UU%tag = trim(wpro%dat_)
    end if

    !! zero the output if add is off
    if(.not. wpro%add) call weidtm_outnull()

    !! setup the adjoint sign
    if(wpro%adj) then
       wpro%opsign = weiB2T
    else
       wpro%opsign = weiT2B
    end if

    if(wpar%nmpi>1) then
       if(wpro%adj) then
          call weiD_mpisprd(UU)
          call weiD_mpizero(DD)
       else
          call weiD_mpisprd(DD)
          call weiD_mpizero(UU)
       end if
    end if

    write(0,'(a,i3,a)') weitab,wpar%impi," begin"
    !===========================================
    if(wpro%outw) then
       weidtm=weidtmof(SLop,WCop,FKop,FXop)
    else 
       weidtm=weidtmoz(SLop,WCop,FKop,FXop)
    end if
    !===========================================
    write(0,'(a,i3,a)') weitab,wpar%impi," end "

    if(wpar%nmpi>1) then
       if(wpro%adj) then
          call weiD_mpigthr(DD)
       else
          call weiD_mpigthr(UU)
       end if
    end if

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
    weidtm=weiOK
  end function weidtm
  !----------------------------------------------------------
  integer function weidtmof(SLop,WCop,FKop,FXop,op5,op6,op7,op8,op9)
    integer,external,optional :: WCop,FKop,FXop,op5,op6,op7,op8,op9
    integer :: izb,iwb, izs,iws,ith
    integer :: st
    integer :: nzbot
#ifdef DEBUG
    character(len=128), parameter :: unit='weidtmof'
#endif
    interface !!----------------------------------------------
       integer function SLop(X) result(st)
         use wei_type
         type(weslow), intent(inout) :: X
       end function SLop
    end interface !!------------------------------------------
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    if(wpro%adj) then
       call weiD_copy(UU,cW)
       wbA: do iwb=1,wpro%nwb
          wpro%bw%o = wpro%nw%o + (iwb-1)*(wpro%bw%n*wpro%bw%d)
          call weiD_read(cW,iwb)
          zbA: do izb=wpro%nzb,1,-1    
             if(wpro%verb) call weipro_followup(izb,iwb)
             call weiS_read(SS,izb)
             st = SLop(SS)

             if(izb==wpro%nzb) then
                nzbot=wpro%nzs-1
             else
                nzbot=wpro%nzs
             end if

             !$OMP  PARALLEL
             !$OMP  DO &
             !$OMP& SCHEDULE(DYNAMIC) &
             !$OMP& PRIVATE(iws,izs,ith)
             do iws=1,wpro%nws
                ith = omp_get_thread_num()+1

                do izs=nzbot,1,-1
                   call weikxmig(cW%bin,iws,izs,ith,WCop,FKop,FXop)
                end do

             end do
             !$OMP END DO
             !$OMP END PARALLEL

          end do zbA
          call weiD_write(cW,iwb)
       end do wbA
       call weiD_add(cW,DD)
    else
       call weiD_copy(DD,cW)
       wbF: do iwb=1,wpro%nwb
          wpro%bw%o = wpro%nw%o + (iwb-1)*(wpro%bw%n*wpro%bw%d)
          call weiD_read(cW,iwb)
          zbF: do izb=1,wpro%nzb       
             if(wpro%verb) call weipro_followup(izb,iwb)
             call weiS_read(SS,izb)
             st = SLop(SS)

             if(izb==wpro%nzb) then
                nzbot=wpro%nzs-1
             else
                nzbot=wpro%nzs
             end if

             !$OMP  PARALLEL
             !$OMP  DO &
             !$OMP& SCHEDULE(DYNAMIC) &
             !$OMP& PRIVATE(iws,izs,ith)
             do iws=1,wpro%nws
                ith = omp_get_thread_num()+1

                do izs=1,nzbot
                   call weikxmig(cW%bin,iws,izs,ith,WCop,FKop,FXop)
                end do

             end do
             !$OMP END DO
             !$OMP END PARALLEL

          end do zbF
          call weiD_write(cW,iwb)
       end do wbF
       call weiD_add(cW,UU)
    end if

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
    weidtmof=weiOK
  end function weidtmof
  !----------------------------------------------------------
  integer function weidtmoz(SLop,WCop,FKop,FXop,op5,op6,op7,op8,op9)
    integer,external,optional :: WCop,FKop,FXop,op5,op6,op7,op8,op9
    integer :: izb,iwb, izs,iws,ith
    integer :: st
    integer :: nzbot
#ifdef DEBUG
    character(len=128), parameter :: unit='weidtmoz'
#endif
    interface !!----------------------------------------------
       integer function SLop(X) result(st)
         use wei_type
         type(weslow), intent(inout) :: X
       end function SLop
    end interface !!------------------------------------------
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    if(wpro%adj) then
       call weiD_copy(UU,cW)
       do izb=wpro%nzb,1,-1
          call weiS_read(SS,izb)
          st = SLop(SS)
          do iwb=1,wpro%nwb
             wpro%bw%o = wpro%nw%o + (iwb-1)*(wpro%bw%n*wpro%bw%d)
             if(wpro%verb) call weipro_followup(izb,iwb)
             call weiD_read(cW,iwb)

             if(izb==wpro%nzb) then
                nzbot=wpro%nzs-1
             else
                nzbot=wpro%nzs
             end if

             !$OMP  PARALLEL
             !$OMP  DO &
             !$OMP& SCHEDULE(DYNAMIC) &
             !$OMP& PRIVATE(iws,izs,ith)
             do iws=1,wpro%nws
                ith = omp_get_thread_num()+1

                do izs=nzbot,1,-1
                   call weikxmig(cW%bin,iws,izs,ith,WCop,FKop,FXop)
                end do

             end do
             !$OMP END DO
             !$OMP END PARALLEL

             call weiD_write(cW,iwb)
          end do
       end do
       call weiD_add(cW,DD)
    else
       call weiD_copy(DD,cW)
       do izb=1,wpro%nzb        
          call weiS_read(SS,izb)
          st = SLop(SS)
          do iwb=1,wpro%nwb
             wpro%bw%o = wpro%nw%o + (iwb-1)*(wpro%bw%n*wpro%bw%d)
             if(wpro%verb) call weipro_followup(izb,iwb)
             call weiD_read(cW,iwb)

             if(izb==wpro%nzb) then
                nzbot=wpro%nzs-1
             else
                nzbot=wpro%nzs
             end if

             !$OMP  PARALLEL
             !$OMP  DO &
             !$OMP& SCHEDULE(DYNAMIC) &
             !$OMP& PRIVATE(iws,izs,ith)
             do iws=1,wpro%nws
                ith = omp_get_thread_num()+1

                do izs=1,nzbot
                   call weikxmig(cW%bin,iws,izs,ith,WCop,FKop,FXop)
                end do

             end do
             !$OMP END DO
             !$OMP END PARALLEL

             call weiD_write(cW,iwb)
          end do
       end do
       call weiD_add(cW,UU)
    end if

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
    weidtmoz=weiOK
  end function weidtmoz
  !---------------------------------------------------------------
end module wei_dtm









