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

contains
  !----------------------------------------------------------------
  subroutine weiW_init(X,data)
    type(wewfld), intent(inout) :: X
    complex, dimension(:,:,:,:,:), pointer :: data

    select case(wpro%operator(1:5))
    case('.Rw.D')
       if(wpar%impi==0) call weiW_par_init(X,data)
       if(wpar%impi/=0) call weiW_scr_init(X,data)
    case('dSwdR')
       if(wpar%impi==0) then
          if(isfile(X%tag)) then
             if(wpro%wTRA) then
                call weiW_his_init(X)
             else
                call weiW_loc_init(X)
             end if
          else
             call seperr('W /= file')
          end if
       end if

       if(wpar%impi/=0) then
          if(wpro%wTRA) then
             call weiW_scr_init(X,data)
          else
             call weiW_loc_init(X)
          end if
       end if
    end select

  end subroutine weiW_init
  !----------------------------------------------------------------
  subroutine weiW_his_init(X)
    type(wewfld), intent(inout) :: X
    character(len=128), parameter :: unit='weiW_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 weiW_report(X)
          call seperr('W.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 weiW_report(X)
          call seperr('W.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 weiW_report(X)
          call seperr('W.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 weiW_report(X)
          call seperr('W.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 weiW_report(X)
          call seperr('W.aw /= pro.aw')
       end if
    end if

    call hinaxis(X%tag,X%az ,6); X%az%l ='z_'
    if(wpro%az%n ==0) then
       wpro%az = X%az
    else
       if(.not. coax( X%az, wpro%az)) then
          call weiW_report(X)
          call seperr('W.az /= pro.az')
       end if
    end if

    call weiW_blockcheck(X)
#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiW_his_init
  !----------------------------------------------------------------
  subroutine weiW_par_init(X,data)
    type(wewfld), intent(inout) :: X
    complex, dimension(:,:,:,:,:), pointer :: data
    character(len=128), parameter :: unit='weiW_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)
    if(wpro%wTRA) then
       X%aw  = wpro%aw;  call houaxis(X%tag,X%aw, 5)
    else !! dont transfer W - allocate just a block
       X%aw  = wpro%nw;  call houaxis(X%tag,X%aw, 5)
    end if
    X%az  = wpro%az;  call houaxis(X%tag,X%az, 6)

    call to_history("esize",8,X%tag)

    call weiW_blockcheck(X)
    call weiW_assign(X,data,cmplx(0.))
#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiW_par_init
  !----------------------------------------------------------------
  subroutine weiW_scr_init(X,data)
    type(wewfld), intent(inout) :: X
    complex, dimension(:,:,:,:,:), pointer :: data
    character(len=128), parameter :: unit='weiW_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)
    X%az  = wpro%az;  call aouaxis(X%tag,X%az, 6)

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

    call weiW_blockcheck(X)
    call weiW_assign(X,data,cmplx(0.))
#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiW_scr_init
  !----------------------------------------------------------------
  subroutine weiW_loc_init(X) !! open local file (dont check much)
    type(wewfld), intent(inout) :: X
    character(len=128), parameter :: unit='weiW_loc_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'
    call hinaxis(X%tag,X%amy,2); X%amy%l='my'
    call hinaxis(X%tag,X%ahx,3); X%ahx%l='hx'
    call hinaxis(X%tag,X%ahy,4); X%ahy%l='hy'
    call hinaxis(X%tag,X%aw ,5); X%aw%l ='w_'
    call hinaxis(X%tag,X%az ,6); X%az%l ='z_'
    
    call weiW_blockcheck(X)

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

    !! DEPTH blocking
    if( (1.0*X%az%n)/wpro%nzs /= X%az%n/wpro%nzs ) then
       call weiW_report(X)
       call seperr('W.az.n is not a multiple of W.nzs')
    end if

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

    X%bites = 8* product( (/wpro%amx%n,wpro%amy%n,wpro%ahx%n,wpro%ahy%n/) )

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiW_blockcheck
  !----------------------------------------------------------------
  subroutine weiW_assign(X,data,v)
    type(wewfld), intent(inout) :: X
    complex, dimension(:,:,:,:,:), pointer :: data
    complex,      intent(in)    :: v
    integer :: izb, izs
    integer :: iwb, iws
    integer :: jmpi
    character(len=128), parameter :: unit='weiW_assign'
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    data=v
    if(wpro%wTRA) then
       if(wpar%impi /= 0) then
          do izb=1,wpro%nzb
             do izs=1,wpro%nzs
                do iwb=1,wpro%nwb
                   do iws=1,wpro%nws
                      call weiW_write(X,data,izb,izs,iwb,iws)
                   end do
                end do
             end do
          end do
       else
          do izb=1,wpro%nzb
             do izs=1,wpro%nzs
                do iwb=1,wpro%nwb * wpar%nmpi
                   do iws=1,wpro%nws
                      call weiW_write(X,data,izb,izs,iwb,iws)
                   end do
                end do
             end do
          end do
       end if
    else
       do izb=1,wpro%nzb
          do izs=1,wpro%nzs
             do iwb=1,wpro%nwb
                do iws=1,wpro%nws
                   call weiW_write(X,data,izb,izs,iwb,iws)
                end do
             end do
          end do
       end do
    end if


#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiW_assign
  !----------------------------------------------------------------  
  subroutine weiW_report(X)
    type(wewfld), intent(in) :: X

    write(0,*)'        ************************************************'
    write(0,*)'        *       ',trim(X%tag),' from ',trim(X%from)
    write(0,'(a,i3,f10.2,f10.3)       ') "         *    Mx: ",X%amx%n, X%amx%o, X%amx%d
    write(0,'(a,i3,f10.2,f10.3)       ') "         *    My: ",X%amy%n, X%amy%o, X%amy%d
    write(0,'(a,i3,f10.2,f10.3)       ') "         *    Hx: ",X%ahx%n, X%ahx%o, X%ahx%d
    write(0,'(a,i3,f10.2,f10.3)       ') "         *    Hy: ",X%ahy%n, X%ahy%o, X%ahy%d
    write(0,'(a,i3,f10.2,f10.3,a,i3,a)') "         *     W: ",X%aw%n,  X%aw%o,  X%aw%d," ->",wpro%nwb," blocks"
    write(0,'(a,i3,f10.2,f10.3,a,i3,a)') "         *     Z: ",X%az%n,  X%az%o,  X%az%d," ->",wpro%nzb," blocks"
    write(0,*)'        ************************************************'
  end subroutine weiW_report
  !--------------------------------------------------------------
  subroutine weiW_read(X,data,izb,izs,iwb,iws)
    type(wewfld), intent(inout) :: X
    complex, dimension(:,:,:,:,:), pointer :: data
    integer, intent(in) :: izb,izs,iwb,iws
    integer :: k
    integer :: stat
    character(len=128), parameter :: unit='weiW_read'
#ifdef DEBUG
    if(wpro%dbg>1) call in(unit)
#endif

    k=(izb-1)*wpro%nzs*X%aw%n &
    + (izs-1)*         X%aw%n &
    + (iwb-1)*wpro%nws &
    +              iws

    stat=sseek_block(X%tag,k-1,X%bites,0)
    call sreed(X%tag, data(:,:,:,:,iws), X%bites)

#ifdef DEBUG
    if(wpro%dbg>1) call out(unit)
#endif
  end subroutine weiW_read
  !----------------------------------------------------------------
  subroutine weiW_write(X,data,izb,izs,iwb,iws)
    type(wewfld), intent(inout) :: X
    complex, dimension(:,:,:,:,:), pointer :: data
    integer, intent(in) :: izb,izs,iwb,iws
    integer :: k
    integer :: stat
    character(len=128), parameter :: unit='weiW_write'
#ifdef DEBUG
    if(wpro%dbg>1) call in(unit)
#endif

    k=(izb-1)*wpro%nzs*X%aw%n &
    + (izs-1)*         X%aw%n &
    + (iwb-1)*wpro%nws &
    +              iws

    stat=sseek_block(X%tag,k-1,X%bites,0)
    call srite(X%tag, data(:,:,:,:,iws), X%bites)

#ifdef DEBUG
    if(wpro%dbg>1) call out(unit)
#endif
  end subroutine weiW_write
  !----------------------------------------------------------------
  subroutine weiW_mpisprd(X,data)
    type(wewfld), intent(inout) :: X
    complex, dimension(:,:,:,:,:), pointer :: data
    integer :: izb, izs
    integer :: iwb, iws

    integer :: mpierr
    integer :: msglen ! message length
    integer :: msgtag ! message tag (arbitrary)
    integer :: msgsou ! message source
    integer :: msgrec ! message receiver
    integer :: jmpi
    character(len=128), parameter :: unit='weiW_mpisprd' 
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    msgtag = 1001
    msglen = size(data)

    if(wpar%impi==0) then
       do jmpi = 1, wpar%nmpi-1
          if(wpro%verb) write(0,'(a,a,i3,a,i3)') weitab,"W send:",wpar%impi,"=>",jmpi
          msgrec = jmpi
          do izb=1,wpro%nzb
             do izs=1,wpro%nzs
                do iwb=1,wpro%nwb

                   do iws=1,wpro%nws
                      call weiW_read (X,data,izb,izs,jmpi*wpro%nwb+iwb,iws)
                   end do
!!$                   call weiW_blkin (X,data,izb,izs,jmpi*wpro%nwb+iwb)
                   CALL MPI_SEND(data,msglen,MPI_COMPLEX,msgrec,msgtag,MPI_COMM_WORLD,       mpierr)
                   
                end do
             end do
          end do
          
       end do
    else
       
       msgsou = 0       
       do izb=1,wpro%nzb
          do izs=1,wpro%nzs
             do iwb=1,wpro%nwb

                CALL MPI_RECV(data,msglen,MPI_COMPLEX,msgsou,msgtag,MPI_COMM_WORLD,status,mpierr)
!!$                call weiW_blkou(X,data,izb,izs,iwb)
                do iws=1,wpro%nws
                   call weiW_write(X,data,izb,izs,iwb,iws)
                end do

             end do
          end do
       end do

    end if

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiW_mpisprd
  !----------------------------------------------------------------
  subroutine weiW_mpigthr(X,data)
    type(wewfld), intent(inout) :: X
    complex, dimension(:,:,:,:,:), pointer :: data
    integer :: izb, izs
    integer :: iwb, iws

    integer :: mpierr
    integer :: msglen ! message length
    integer :: msgtag ! message tag (arbitrary)
    integer :: msgsou ! message source
    integer :: msgrec ! message receiver
    integer :: jmpi
    complex, pointer :: tmp(:,:,:,:,:)
    character(len=128), parameter :: unit='weiW_mpigthr' 
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    msgtag = 1001
    msglen = size(data)

    if(wpar%impi==0) then
       call weimem_add(8,size(data,1)*size(data,2)*size(data,3)*size(data,4)*size(data,5),unit)
       allocate(     tmp(size(data,1),size(data,2),size(data,3),size(data,4),size(data,5)))

       do jmpi = 1, wpar%nmpi-1
          if(wpro%verb) write(0,'(a,a,i3,a,i3)') weitab,"W recv:",wpar%impi,"<=",jmpi

          msgsou = jmpi
          do izb=1,wpro%nzb
             do izs=1,wpro%nzs
                do iwb=1,wpro%nwb

                   do iws=1,wpro%nws
                      call weiW_read (X,data,izb,izs,jmpi*wpro%nwb+iwb,iws)
                   end do
!!$                   call weiW_blkin(X,data,izb,izs,jmpi*wpro%nwb+iwb)

                   CALL MPI_RECV(  tmp,msglen,MPI_COMPLEX,msgsou,msgtag,MPI_COMM_WORLD,status,mpierr)
                   data = data + tmp

!!$                   call weiW_blkou(X,data,izb,izs,jmpi*wpro%nwb+iwb)
                   do iws=1,wpro%nws
                      call weiW_write(X,data,izb,izs,jmpi*wpro%nwb+iwb,iws)
                   end do

                end do
             end do
          end do

       end do

       call weimem_cut(8,size(data,1)*size(data,2)*size(data,3)*size(data,4)*size(data,5),unit)
       deallocate(tmp)
    else

       msgrec = 0
       do izb=1,wpro%nzb
          do izs=1,wpro%nzs
             do iwb=1,wpro%nwb

                do iws=1,wpro%nws
                   call weiW_read(X,data,izb,izs,iwb,iws)
                end do
!!$                call weiW_blkin(X,data,izb,izs,iwb)
                CALL MPI_SEND( data,msglen,MPI_COMPLEX,msgrec,msgtag,MPI_COMM_WORLD,       mpierr)

             end do
          end do
       end do

    end if

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiW_mpigthr
  !----------------------------------------------------------------
  subroutine weiW_mpizero(X,data)
    type(wewfld), intent(inout) :: X
    complex, dimension(:,:,:,:,:), pointer :: data
    character(len=128), parameter :: unit='weiW_mpizero' 
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif
    if(wpar%impi/=0) call weiW_assign(X,data,cmplx(0.))
#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiW_mpizero
  !----------------------------------------------------------------
  subroutine weiW_blkin(X,data,izb,izs,iwb)
    type(wewfld), intent(inout) :: X
    complex, dimension(:,:,:,:,:), pointer :: data
    integer, intent(in) :: izb,izs,iwb
    integer :: k
    integer :: stat
    character(len=128), parameter :: unit='weiW_blkin'
#ifdef DEBUG
    if(wpro%dbg>1) call in(unit)
#endif

    k=(izb-1)*wpro%nzs*X%aw%n &
    + (izs-1)*         X%aw%n &
    + (iwb-1)*wpro%nws &
    + 1

    stat=sseek_block(X%tag,k-1 ,wpro%nws*X%bites,0)
    call       sreed(X%tag,data,wpro%nws*X%bites)

#ifdef DEBUG
    if(wpro%dbg>1) call out(unit)
#endif
  end subroutine weiW_blkin
  !----------------------------------------------------------------
  subroutine weiW_blkou(X,data,izb,izs,iwb)
    type(wewfld), intent(inout) :: X
    complex, dimension(:,:,:,:,:), pointer :: data
    integer, intent(in) :: izb,izs,iwb
    integer :: k
    integer :: stat
    character(len=128), parameter :: unit='weiW_blkou'
#ifdef DEBUG
    if(wpro%dbg>1) call in(unit)
#endif

    k=(izb-1)*wpro%nzs*X%aw%n &
    + (izs-1)*         X%aw%n &
    + (iwb-1)*wpro%nws &
    + 1
    
    stat=sseek_block(X%tag,k-1 ,wpro%nws*X%bites,0)
    call       srite(X%tag,data,wpro%nws*X%bites)

#ifdef DEBUG
    if(wpro%dbg>1) call out(unit)
#endif
  end subroutine weiW_blkou
  !----------------------------------------------------------------
end module wei_W





