module wei_D
  use wei_util
  use wei_type
  use wei_pro
  use wei_par
  !---------------
  implicit none

contains
  !--------------------------------------------------------------
  subroutine weiD_init(X)
    type(wedata), intent(inout) :: X

    if(wpar%impi==0) then
       if(isfile(X%tag)) then
          call weiD_his_init(X)
       else
          call weiD_par_init(X)
       end if
    end if

    if(wpar%impi/=0) call weiD_scr_init(X)

  end subroutine weiD_init
  !----------------------------------------------------------------
  subroutine weiD_aux_init(X)
    type(wedata),       intent(inout) :: X
    character(len=256) :: buffer

    write(buffer, '(a,a,a)') ".",trim(X%tag), trim(wpar%suffix)
    X%tag = trim(buffer(1:128))
    call weiD_scr_init(X)

  end subroutine weiD_aux_init
  !----------------------------------------------------------------
  subroutine weiD_his_init(X)
    type(wedata), intent(inout) :: X
    character(len=128), parameter :: unit='weiD_his_init'
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    call auxinout(X%tag)
    X%from = "history"

    call hinaxis(X%tag,X%amx,1); X%amx%l='mx'
    if(wpro%amx%n ==0) then
       wpro%amx = X%amx
    else
       if(.not. coax( X%amx, wpro%amx)) then
          call weiD_report(X)
          call seperr('D.amx /= pro.amx')
       end if
    end if

    call hinaxis(X%tag,X%amy,2); X%amy%l='my'
    if(wpro%amy%n ==0) then
       wpro%amy = X%amy
    else
       if(.not. coax( X%amy, wpro%amy)) then
          call weiD_report(X)
          call seperr('D.amy /= pro.amy')
       end if
    end if

    call hinaxis(X%tag,X%ahx,3); X%ahx%l='hx'
    if(wpro%ahx%n ==0) then
       wpro%ahx = X%ahx
    else
       if(.not. coax( X%ahx, wpro%ahx)) then
          call weiD_report(X)
          call seperr('D.ahx /= pro.ahx')
       end if
    end if

    call hinaxis(X%tag,X%ahy,4); X%ahy%l='hy'
    if(wpro%ahy%n ==0) then
       wpro%ahy = X%ahy
    else
       if(.not. coax( X%ahy, wpro%ahy)) then
          call weiD_report(X)
          call seperr('D.ahy /= pro.ahy')
       end if
    end if

    call hinaxis(X%tag,X%aw ,5); X%aw%l ='w_'
    if(wpro%aw%n ==0) then
       wpro%aw = X%aw
    else
       if(.not. coax( X%aw, wpro%aw)) then
          call weiD_report(X)
          call seperr('D.aw /= pro.aw')
       end if
    end if

    call weiD_blockcheck(X)
    call weiD_allocate(X)

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiD_his_init
  !----------------------------------------------------------------
  subroutine weiD_par_init(X)
    type(wedata), intent(inout) :: X
    character(len=128), parameter :: unit='weiD_par_init' 
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

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

    X%amx = wpro%amx; call houaxis(X%tag,X%amx,1)
    X%amy = wpro%amy; call houaxis(X%tag,X%amy,2)
    X%ahx = wpro%ahx; call houaxis(X%tag,X%ahx,3)
    X%ahy = wpro%ahy; call houaxis(X%tag,X%ahy,4)
    X%aw  = wpro%aw ; call houaxis(X%tag,X%aw, 5)

    call auxputch("esize","i",8,X%tag)

    call weiD_blockcheck(X)
    call weiD_allocate(X)
    call weiD_assign(X,cmplx(0.))

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiD_par_init
  !----------------------------------------------------------------
  subroutine weiD_scr_init(X)
    type(wedata), intent(inout) :: X
    character(len=128), parameter :: unit='weiD_scr_init' 
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

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

    X%amx = wpro%amx; call aouaxis(X%tag,X%amx,1)
    X%amy = wpro%amy; call aouaxis(X%tag,X%amy,2)
    X%ahx = wpro%ahx; call aouaxis(X%tag,X%ahx,3)
    X%ahy = wpro%ahy; call aouaxis(X%tag,X%ahy,4)
    X%aw  = wpro%nw;  call aouaxis(X%tag,X%aw, 5)

    call auxputch("esize","i",8,X%tag)

    call weiD_blockcheck(X)
    call weiD_allocate(X)
    call weiD_assign(X,(0.,0.))
#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiD_scr_init
  !----------------------------------------------------------------
  subroutine weiD_blockcheck(X)
    type(wedata), intent(inout) :: X

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

  end subroutine weiD_blockcheck
  !----------------------------------------------------------------
  subroutine weiD_allocate(X)
    type(wedata), intent(inout) :: X
    integer :: er
    character(len=128), parameter :: unit='weiD_allocate'
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

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

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiD_allocate
  !----------------------------------------------------------------
  subroutine weiD_assign(X,v)
    type(wedata), intent(inout) :: X
    complex, intent(in) :: v
    integer :: iwb
    integer :: jmpi
    character(len=128), parameter :: unit='weiD_assign' 
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    X%bin=v
    if(wpar%impi /= 0) then
       do iwb=1,wpro%nwb
          call weiD_write(X,iwb)
       end do
    else
       do jmpi=0,wpar%nmpi-1
          do iwb=1,wpro%nwb
             call weiD_write(X,jmpi*wpro%nwb+iwb)
          end do
       end do
    end if

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiD_assign
  !----------------------------------------------------------------
  subroutine weiD_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%n,")",  X%aw%o,  X%aw%d," ->",wpro%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 weiD_report
  !--------------------------------------------------------------
  subroutine weiD_read(X,k)
    type(wedata), intent(inout) :: X
    integer,      intent(in)    :: k
    integer :: stat
    character(len=128), parameter :: unit='weiD_read'
#ifdef DEBUG
    if(wpro%dbg>1) call in(unit)
#endif

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

#ifdef DEBUG
    if(wpro%dbg>1) call out(unit)
#endif
  end subroutine weiD_read
  !----------------------------------------------------------------
  subroutine weiD_write(X,k)
    type(wedata), intent(inout) :: X
    integer,      intent(in)    :: k
    integer :: stat
    character(len=128), parameter :: unit='weiD_write' 
#ifdef DEBUG
    if(wpro%dbg>1) call in(unit)
#endif

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

#ifdef DEBUG
    if(wpro%dbg>1) call out(unit)
#endif
  end subroutine weiD_write
  !----------------------------------------------------------------
  subroutine weiD_copy(X,Y)
    type(wedata), intent(inout) :: X
    type(wedata), intent(inout) :: Y
    integer :: iwb
    character(len=128), parameter :: unit='weiD_copy'
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    do iwb=1,wpro%nwb
       call weiD_read (X,iwb)
       Y%bin = X%bin
       call weiD_write(Y,iwb)
    end do

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiD_copy
  !----------------------------------------------------------------
  subroutine weiD_add(X,Y)
    type(wedata), intent(inout) :: X
    type(wedata), intent(inout) :: Y
    integer :: iwb
    character(len=128), parameter :: unit='weiD_add'
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    do iwb=1,wpro%nwb
       call weiD_read (X,iwb)
       call weiD_read (Y,iwb)
       Y%bin = Y%bin + X%bin
       call weiD_write(Y,iwb)
    end do

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiD_add
  !----------------------------------------------------------------  
  subroutine weiD_mpisprd(X)
    type(wedata), intent(inout) :: X
    integer :: mpierr
    integer :: msglen ! message length
    integer :: msgtag ! message tag (arbitrary)
    integer :: msgsou ! message source
    integer :: msgrec ! message receiver
    integer :: iwb,jmpi
    character(len=128), parameter :: unit='weiD_mpisprd' 
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    msgtag = 1001
    msglen = size(X%bin,1)*size(X%bin,2)*size(X%bin,3)*size(X%bin,4)*size(X%bin,5)

    if(wpar%impi==0) then
       do jmpi = 1, wpar%nmpi-1
          if(wpro%verb) write(0,'(a,a,i3,a,i3)') weitab,"D send:",wpar%impi,"=>",jmpi
          msgrec = jmpi
          do iwb=1,wpro%nwb
             call weiD_read(X,jmpi*wpro%nwb+iwb)
             CALL MPI_SEND(X%bin,msglen,MPI_COMPLEX,msgrec,msgtag,MPI_COMM_WORLD,       mpierr)
          end do
       end do
    else
       msgsou = 0
       do iwb=1,wpro%nwb
          ;  CALL MPI_RECV(X%bin,msglen,MPI_COMPLEX,msgsou,msgtag,MPI_COMM_WORLD,status,mpierr)
          ;  call weiD_write(X,iwb)
       end do
    end if

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiD_mpisprd
  !----------------------------------------------------------------
  subroutine weiD_mpigthr(X)
    type(wedata), intent(inout) :: X
    integer :: mpierr
    integer :: msglen ! message length
    integer :: msgtag ! message tag (arbitrary)
    integer :: msgsou ! message source
    integer :: msgrec ! message receiver
    integer :: iwb,jmpi
    complex, pointer :: tmp(:,:,:,:,:)
    character(len=128), parameter :: unit='weiD_mpigthr' 
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    msgtag = 1001
    msglen = size(X%bin,1)*size(X%bin,2)*size(X%bin,3)*size(X%bin,4)*size(X%bin,5)
    if(wpar%impi==0) then
       call weimem_add(8,size(X%bin,1)*size(X%bin,2)*size(X%bin,3)*size(X%bin,4)*size(X%bin,5),unit)
       allocate(     tmp(size(X%bin,1),size(X%bin,2),size(X%bin,3),size(X%bin,4),size(X%bin,5)))

       do jmpi = 1, wpar%nmpi-1
          if(wpro%verb) write(0,'(a,a,i3,a,i3)') weitab,"D recv:",wpar%impi,"<=",jmpi
          msgsou = jmpi
          do iwb=1,wpro%nwb
             CALL MPI_RECV(  tmp,msglen,MPI_COMPLEX,msgsou,msgtag,MPI_COMM_WORLD,status,mpierr)

             call weiD_read (X,jmpi*wpro%nwb+iwb)
             X%bin = X%bin + tmp
             call weiD_write(X,jmpi*wpro%nwb+iwb)
          end do

       end do

       call weimem_cut(8,size(X%bin,1)*size(X%bin,2)*size(X%bin,3)*size(X%bin,4)*size(X%bin,5),unit)
       deallocate(tmp)
    else
       msgrec = 0
       do iwb=1,wpro%nwb
          ;  call weiD_read(X,iwb)
          ;  CALL MPI_SEND(X%bin,msglen,MPI_COMPLEX,msgrec,msgtag,MPI_COMM_WORLD,       mpierr)
       end do
    end if

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiD_mpigthr
  !----------------------------------------------------------------
  subroutine weiD_mpizero(X)
    type(wedata), intent(inout) :: X
    character(len=128), parameter :: unit='weiD_mpizero' 
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    if(wpar%impi/=0) call weiD_assign(X,cmplx(0.))

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiD_mpizero
  !----------------------------------------------------------------
end module wei_D











