module wei_pro
  use sep
  use wei_util
  use wei_external
  use wei_par
  !---------------
  implicit none
  !----------------------------------------------------------------
  type weipro
     character(len=128) :: operator

     logical            :: adj, add
     integer            :: opsign
     !! simultate conjg(phase) to avoid 'if'
     !! opsign=weiB2T (e.g.  modeling)
     !! opsign=weiT2B (e.g. migration)

     logical            :: verb    !! verbose
     integer            :: dbg     !! debug flag

     logical            :: outw    !! frequency in outer loop
     logical            :: pCIG    !! ph images
     logical            :: wTRA    !! wavefield transfer

     character(len=128) :: mod_
     character(len=128) :: dat_

     integer            :: nzs     !! number of depth     slices in a block
     integer            :: nws     !! number of frequency slices in a block

     integer            :: nzb     !! number of depth     blocks
     integer            :: nwb     !! number of frequency blocks

     integer            :: nfk     !! number of f-k ops

     real               :: pi

     type(axis)         :: amx,amy !! midpoint
     type(axis)         :: ahx,ahy !! offset
     type(axis)         :: apx,apy !! offset ray parameter
     type(axis)         :: az,bz   !! depth

     type(axis)         :: aw      !! all   frequencies
     type(axis)         :: nw      !! node  frequencies
     type(axis)         :: bw      !! block frequencies

     logical            :: memverb !! print memory allocation info
     real               :: mem     !! allocated memory (in Mb)
  end type weipro
  !----------------------------------------------------------------

  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!
  type(weipro), public :: wpro
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!

  integer, public :: weiOK=1
  integer, public :: weiB2T=+1,weiT2B=-1
  !! go Top2Bottom or Bottom2Top

contains
  !----------------------------------------------------------------
  subroutine weipro_init()
    character(len=128) :: operator =""
    character(len=128) :: operation=""
    logical            :: adj,add
    integer            :: ierr
#ifdef DEBUG
    character(len=128), parameter :: unit='weipro_init'
#endif

    call from_param('dbg',  wpro%dbg,       0 )
    call from_param('verb', wpro%verb, .false.)
    if(wpro%dbg>0) wpro%verb =.true.

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

    !! processing flags
    call from_param('adj', wpro%adj, .false.)
    call from_param('add', wpro%add, .false.)
    call from_param("outw",wpro%outw,.true. )
    call from_param("pCIG",wpro%pCIG,.false.)
    call from_param("wTRA",wpro%wTRA,.true. )

    call from_param('operator',operator,'.....')
    wpro%operator = operator !! .R_.D;.Rw.D;dS_dR;dSwdR

    call from_param('operation',operation,'dottest...')
    if(operation(1:10) .eq. 'inversion.') wpro%adj=.true. 

    !! number of slices in a block slices
    call from_param("nzs",wpro%nzs,1)
    call from_param("nws",wpro%nws,1); wpro%nwb=0;

    !! axes info from the parameter file
    !! amx,amy,ahx,ahy,apx,apy, aw,az, bw,bz
    !! here nws and nzs must be known
    call weipro_parameters()

    wpro%nzb = wpro%az%n             / wpro%nzs
    wpro%nwb = wpro%aw%n / wpar%nmpi / wpro%nws

    !! no of reference velocities
    call from_param("nfk",wpro%nfk,1)

    wpro%pi=acos(-1.)

    !! adjust w block axes for the various nodes
    !! frequencies allocated to each node
    wpro%nw%n = wpro%aw%n / wpar%nmpi
    wpro%nw%d = wpro%aw%d
    wpro%nw%o = wpro%aw%o + wpar%impi * (wpro%aw%n/wpar%nmpi) * wpro%aw%d

    call from_param("memverb",wpro%memverb,.false.)
    wpro%mem = 0.

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weipro_init
  !----------------------------------------------------------------
  subroutine weipro_parameters()

    wpro%amx%l='mx'; call pinaxis(wpro%amx)
    wpro%amy%l='my'; call pinaxis(wpro%amy)

    wpro%ahx%l='hx'; call pinaxis(wpro%ahx)
    wpro%ahy%l='hy'; call pinaxis(wpro%ahy)

    wpro%apx%l='px'; call pinaxis(wpro%apx)
    wpro%apy%l='py'; call pinaxis(wpro%apy)

    wpro%aw%l ='w_'; call pinaxis(wpro%aw )
    wpro%az%l ='z_'; call pinaxis(wpro%az )

    wpro%bz = wpro%az; wpro%bz%n=wpro%nzs
    wpro%bw = wpro%aw; wpro%bw%n=wpro%nws

  end subroutine weipro_parameters
  !----------------------------------------------------------------
  subroutine weipro_bcast()
    integer :: mpierr

    call MPI_BCAST(wpro%operator,128,MPI_CHARACTER,0,MPI_COMM_WORLD,mpierr)

    call MPI_BCAST(wpro%adj, 1,MPI_LOGICAL,0,MPI_COMM_WORLD,mpierr) ! adj
    call MPI_BCAST(wpro%add, 1,MPI_LOGICAL,0,MPI_COMM_WORLD,mpierr) ! add
    call MPI_BCAST(wpro%verb,1,MPI_LOGICAL,0,MPI_COMM_WORLD,mpierr) ! verb

    call MPI_BCAST(wpro%amx%n,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) ! amx
    call MPI_BCAST(wpro%amx%o,1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! amx
    call MPI_BCAST(wpro%amx%d,1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! amx

    call MPI_BCAST(wpro%amy%n,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) ! amy
    call MPI_BCAST(wpro%amy%o,1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! amy
    call MPI_BCAST(wpro%amy%d,1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! amy

    call MPI_BCAST(wpro%ahx%n,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) ! ahx
    call MPI_BCAST(wpro%ahx%o,1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! ahx
    call MPI_BCAST(wpro%ahx%d,1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! ahx

    call MPI_BCAST(wpro%ahy%n,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) ! ahy
    call MPI_BCAST(wpro%ahy%o,1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! ahy
    call MPI_BCAST(wpro%ahy%d,1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! ahy

    call MPI_BCAST(wpro%apx%n,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) ! apx
    call MPI_BCAST(wpro%apx%o,1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! apx
    call MPI_BCAST(wpro%apx%d,1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! apx

    call MPI_BCAST(wpro%apy%n,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) ! apy
    call MPI_BCAST(wpro%apy%o,1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! apy
    call MPI_BCAST(wpro%apy%d,1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! apy

    call MPI_BCAST(wpro%az%n, 1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) ! az
    call MPI_BCAST(wpro%az%o, 1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! az
    call MPI_BCAST(wpro%az%d, 1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! az

    call MPI_BCAST(wpro%aw%n, 1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) ! aw
    call MPI_BCAST(wpro%aw%o, 1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! aw
    call MPI_BCAST(wpro%aw%d, 1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! aw

    call MPI_BCAST(wpro%nw%n, 1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) ! nw
    call MPI_BCAST(wpro%nw%o, 1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! nw
    call MPI_BCAST(wpro%nw%d, 1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! nw

    call MPI_BCAST(wpro%bw%n, 1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) ! bw
    call MPI_BCAST(wpro%bw%o, 1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! bw
    call MPI_BCAST(wpro%bw%d, 1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! bw

    call MPI_BCAST(wpro%outw, 1,MPI_LOGICAL,0,MPI_COMM_WORLD,mpierr) ! outw
    call MPI_BCAST(wpro%pCIG, 1,MPI_LOGICAL,0,MPI_COMM_WORLD,mpierr) ! pCIG
    call MPI_BCAST(wpro%wTRA, 1,MPI_LOGICAL,0,MPI_COMM_WORLD,mpierr) ! wTRA

    call MPI_BCAST(wpro%nfk,  1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) ! nfk
    call MPI_BCAST(wpro%nzs,  1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) ! nzs
    call MPI_BCAST(wpro%nws,  1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) ! nws

    call MPI_BCAST(wpro%nzb,  1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) ! nzb
    call MPI_BCAST(wpro%nwb,  1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) ! nwb

    call MPI_BCAST(wpro%pi ,  1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr) ! pi

    !! adjust w block axes for the various nodes
    !! frequencies allocated to each node
    wpro%nw%n = wpro%aw%n / wpar%nmpi
    wpro%nw%d = wpro%aw%d
    wpro%nw%o = wpro%aw%o + wpar%impi * (wpro%aw%n/wpar%nmpi) * wpro%aw%d

    call MPI_BCAST(wpro%memverb,1,MPI_LOGICAL,0,MPI_COMM_WORLD,mpierr) ! memverb

  end subroutine weipro_bcast
  !----------------------------------------------------------------
  subroutine weipro_reset(adj,add,mod_,dat_)
    logical,            intent(in) :: adj,add
    character(len=128), intent(in) :: mod_,dat_

    wpro%adj  = adj
    wpro%add  = add
    wpro%mod_ = mod_
    wpro%dat_ = dat_
  end subroutine weipro_reset
  !----------------------------------------------------------------  
  subroutine weipro_report()
    call weiline()
    write(0,'(a,a,i2)') weitab," WEIPRO REPORT", wpar%impi
    call weiline()
    write(0,'(a,a,a)')  weitab," operator=",wpro%operator (1: 5)
    call weiline()
    !    write(0,'(a,a,l1,a,l1,a,l1,a,i2)') weitab," adj=",wpro%adj," add=",wpro%add &
    !    ,                 " verb=",wpro%verb," dbg=",wpro%dbg

    write(0,'(a,a,l1,a,l1)')      weitab," adj=",wpro%adj,"  add=",wpro%add
    write(0,'(a,a,l1,a,i2,a,l1)') weitab,"verb=",wpro%verb," dbg=",wpro%dbg,"memverb=",wpro%memverb

    call weiline()
    write(0,'(a,a,i4,f13.5,f13.5)') weitab," z = ", wpro%az%n ,  wpro%az%o ,  wpro%az%d
    write(0,'(a,a,i4,f13.5,f13.5)') weitab," w = ", wpro%aw%n ,  wpro%aw%o ,  wpro%aw%d
    write(0,'(a,a,i4,f13.5,f13.5)') weitab," w = ", wpro%nw%n ,  wpro%nw%o ,  wpro%nw%d
    write(0,'(a,a,i4,f13.5,f13.5)') weitab," mx= ", wpro%amx%n,  wpro%amx%o,  wpro%amx%d
    write(0,'(a,a,i4,f13.5,f13.5)') weitab," my= ", wpro%amy%n,  wpro%amy%o,  wpro%amy%d
    write(0,'(a,a,i4,f13.5,f13.5)') weitab," hx= ", wpro%ahx%n,  wpro%ahx%o,  wpro%ahx%d
    write(0,'(a,a,i4,f13.5,f13.5)') weitab," hy= ", wpro%ahy%n,  wpro%ahy%o,  wpro%ahy%d
    write(0,'(a,a,i4,f13.5,f13.5)') weitab," px= ", wpro%apx%n,  wpro%apx%o,  wpro%apx%d
    write(0,'(a,a,i4,f13.5,f13.5)') weitab," py= ", wpro%apy%n,  wpro%apy%o,  wpro%apy%d
    call weiline()
    write(0,'(a,a,l1,a,l1,a,l1)')  weitab," outw=",wpro%outw," pCIG=",wpro%pCIG," wTRA=", wpro%wTRA
    call weiline()
    write(0,'(a,a,i3,a,i3)') weitab," nzs=",wpro%nzs," nzb=",wpro%nzb
    write(0,'(a,a,i3,a,i3)') weitab," nws=",wpro%nws," nwb=",wpro%nwb
    write(0,'(a,a,i3     )') weitab," nfk=",wpro%nfk

!!$    write(0,'(a,a,i3,a,i3,a,i3,a,i3,a,i2)') weitab," nzs=",wpro%nzs," nws=",wpro%nws&
!!$    ,                                              " nzb=",wpro%nzb," nwb=",wpro%nwb&
!!$    ,                                              " nfk=",wpro%nfk

    call weiline()
    write(0,'(a,a,a,a,a)')  weitab," mod=",trim(wpro%mod_)," dat=",trim(wpro%dat_)
    call weiline()

  end subroutine weipro_report
  !----------------------------------------------------------------
  subroutine weipro_follow()

    call weiline()
    if(wpro%adj) then
       if(wpro%add) then
          write(0,'(a,a,a,a,a,a,a)') weitab,'[',trim(wpro%operator),']  ',trim(wpro%mod_),'<+<',trim(wpro%dat_)
       else
          write(0,'(a,a,a,a,a,a,a)') weitab,'[',trim(wpro%operator),']  ',trim(wpro%mod_),'< <',trim(wpro%dat_)
       end if
    else
       if(wpro%add) then
          write(0,'(a,a,a,a,a,a,a)') weitab,'[',trim(wpro%operator),']  ',trim(wpro%mod_),'>+>',trim(wpro%dat_)
       else
          write(0,'(a,a,a,a,a,a,a)') weitab,'[',trim(wpro%operator),']  ',trim(wpro%mod_),'> >',trim(wpro%dat_)
       end if
    end if
    call weiline()
  end subroutine weipro_follow
  !----------------------------------------------------------------
  subroutine weipro_followup(izb,iwb)
    integer, intent(in)    ::izb,iwb
    integer :: i
    integer :: allblocks, doneblocks
    real    :: completed

    allblocks = wpro%nzb*wpro%nwb

    select case(wpro%operator(1:5))
    case('.R_.D','.Rw.D')

       if(wpro%adj) then !! migration
          if(wpro%outw) then
             doneblocks= (iwb-1)*wpro%nzb + izb -1
          else
             doneblocks= (izb-1)*wpro%nwb + iwb -1
          end if
       else             !! modeling
          if(wpro%outw) then
             doneblocks= (iwb-1)*wpro%nzb + wpro%nzb-izb
          else
             doneblocks= (wpro%nzb-izb)*wpro%nwb + iwb -1
          end if
       end if

    case('dS_dR','dSwdR','.D_.U')

       if(wpro%adj) then !! backprojection
          if(wpro%outw) then
             doneblocks= (iwb-1)*wpro%nzb + wpro%nzb-izb
          else
             doneblocks= (wpro%nzb-izb)*wpro%nwb + iwb -1
          end if
       else             !! scattering
          if(wpro%outw) then 
             doneblocks= (iwb-1)*wpro%nzb + izb -1
          else
             doneblocks= (izb-1)*wpro%nwb + iwb -1
          end if
       end if

    end select

    completed = 100.0 * doneblocks/allblocks !! percent completed

    !    write(0,'(a,a,f5.2,a,i4,a,i4,a,i4,a,i4)') weitab,"done "&
    !    , completed," %       (z)",izb,"/",wpro%nzb &
    !    ,             "       (w)",iwb,"/",wpro%nwb

    write(0,'(a,i3,a,f5.2,a,i4,a,i4,a,i4,a,i4)') &
    weitab,wpar%impi," done ", completed," %     " &
    , izb,"  <z>",wpro%nzb,"      " &
    , iwb,"  <w>",wpro%nwb

  end subroutine weipro_followup
  !----------------------------------------------------------------
  subroutine weipro_radians()
#ifdef DEBUG
    character(len=128), parameter :: unit='weipro_radians'
    if(wpro%dbg>0) call in(unit)
#endif
    wpro%aw%o = 2*wpro%pi * wpro%aw%o
    wpro%aw%d = 2*wpro%pi * wpro%aw%d

    wpro%nw%o = 2*wpro%pi * wpro%nw%o
    wpro%nw%d = 2*wpro%pi * wpro%nw%d

    wpro%bw%o = 2*wpro%pi * wpro%bw%o
    wpro%bw%d = 2*wpro%pi * wpro%bw%d
#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weipro_radians
  !----------------------------------------------------------------
  subroutine weipro_check_hp()
#ifdef DEBUG
    character(len=128), parameter :: unit='weipro_check_hp'
    if(wpro%dbg>0) call in(unit)
#endif
    if(wpro%apx%n /= wpro%ahx%n) then
       call seperr('apx.n /= ahx.n')
    end if
    if(wpro%apy%n /= wpro%ahy%n) then
       call seperr('apy.n /= ahy.n')
    end if
#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weipro_check_hp
  !----------------------------------------------------------------
  subroutine weimem(caller)
    character(len=*), intent(in), optional :: caller
    if(present(caller)) then
       write(0,'(a,i3,a,f10.4,a,a,a)') weitab, wpar%impi," mem= ",wpro%mem,' Mb', '  @ ',trim(caller)
    else
       write(0,'(a,i3,a,f10.4,a    )') weitab, wpar%impi," mem= ",wpro%mem,' Mb'
    end if
  end subroutine weimem
  !----------------------------------------------------------------
  subroutine weimem_add(   esize,elems,caller)
    integer, intent(in) :: esize,elems
    character(len=*), intent(in), optional :: caller
    real                :: scale
    scale = esize/1024./1024. !! convert to Mb
    wpro%mem = wpro%mem + scale * elems

    if(wpro%memverb) call weimem(caller)
  end subroutine weimem_add
  !----------------------------------------------------------------
  subroutine weimem_cut(   esize,elems,caller)
    integer, intent(in) :: esize,elems
    character(len=*), intent(in), optional :: caller
    real                :: scale
    scale = esize/1024./1024. !! convert to Mb
    wpro%mem = wpro%mem - scale * elems

    if(wpro%memverb) call weimem(caller)
  end subroutine weimem_cut
  !----------------------------------------------------------------
end module wei_pro






































