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

contains
  !----------------------------------------------------------------
  subroutine weiR_init(X)
    type(weimag), intent(inout) :: X

    if(wpar%impi==0) then
       if(isfile(X%tag)) then
          call weiR_his_init(X)
       else
          call weiR_par_init(X)
       end if
    end if

    if(wpar%impi/=0) call weiR_scr_init(X)
  end subroutine weiR_init
  !----------------------------------------------------------------
  subroutine weiR_his_init(X)
    type(weimag), intent(inout) :: X
    character(len=128), parameter :: unit='weiR_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 weiR_report(X)
          call seperr('R.amx /= wpro.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 weiR_report(X)
          call seperr('R.amy /= wpro.amy')
       end if
    end if

    if(wpro%pCIG) then
       call hinaxis(X%tag,X%apx,3); X%apx%l='px'
       if(wpro%apx%n ==0) then
          wpro%apx = X%apx
       else
          if(.not. coax( X%apx, wpro%apx)) then
             call weiR_report(X)
             call seperr('R.apx /= wpro.apx')
          end if
       end if

       call hinaxis(X%tag,X%apy,4); X%apy%l='py'
       if(wpro%apy%n ==0) then
          wpro%apy = X%apy
       else
          if(.not. coax( X%apy, wpro%apy)) then
             call weiR_report(X)
             call seperr('R.apy /= wpro.apy')
          end if
       end if
    else
       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 weiR_report(X)
             call seperr('R.ahx /= wpro.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 weiR_report(X)
             call seperr('R.ahy /= wpro.ahy')
          end if
       end if
    end if

    call hinaxis(X%tag,X%az ,5); X%az%l ='z_'
    if(wpro%az%n ==0) then
       wpro%az = X%az
    else
       if(.not. coax( X%az, wpro%az)) then
          call weiR_report(X)
          call seperr('R.az /= wpro.az')
       end if
    end if

    call weiR_blockcheck(X)
    call weiR_allocate(X)

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

    call auxscr(trim(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)
    if(wpro%pCIG) then
       X%apx = wpro%apx; call houaxis(X%tag,X%apx,3)
       X%apy = wpro%apy; call houaxis(X%tag,X%apy,4)
    else
       X%ahx = wpro%ahx; call houaxis(X%tag,X%ahx,3)
       X%ahy = wpro%ahy; call houaxis(X%tag,X%ahy,4)
    end if
    X%az  = wpro%az;     call houaxis(X%tag,X%az, 5)

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

    call weiR_blockcheck(X)
    call weiR_allocate(X)
    call weiR_assign(X,cmplx(0.))

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

    call auxscr(trim(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)
    if(wpro%pCIG) then
       X%apx = wpro%apx; call aouaxis(X%tag,X%apx,3)
       X%apy = wpro%apy; call aouaxis(X%tag,X%apy,4)
    else
       X%ahx = wpro%ahx; call aouaxis(X%tag,X%ahx,3)
       X%ahy = wpro%ahy; call aouaxis(X%tag,X%ahy,4)
    end if
    X%az  = wpro%az;     call aouaxis(X%tag,X%az, 5)

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

    call weiR_blockcheck(X)
    call weiR_allocate(X)
    call weiR_assign(X,cmplx(0.))

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiR_scr_init
  !----------------------------------------------------------------
  subroutine weiR_blockcheck(X)
    type(weimag), intent(inout) :: X

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

  end subroutine weiR_blockcheck
  !--------------------------------------------------------------
  subroutine weiR_allocate(X)
    type(weimag), intent(inout) :: X
    integer :: er
    character(len=128), parameter :: unit='weiR_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%nzs,unit)
    !!--------------------
    allocate(X%bin(wpro%amx%n,wpro%amy%n,wpro%ahx%n,wpro%ahy%n,wpro%nzs),stat=er)
    if(er/=0) call seperr('cant allocate R.bin')
    X%bin=(0.0,0.0)
    !!--------------------
    X%bites= 8* product((/wpro%amx%n,wpro%amy%n,wpro%ahx%n,wpro%ahy%n,wpro%nzs/))

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiR_allocate
  !----------------------------------------------------------------
  subroutine weiR_assign(X,v)
    type(weimag), intent(inout) :: X
    complex,      intent(in)    :: v
    integer :: izb
    character(len=128), parameter :: unit='weiR_assign'
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    X%bin=v
    do izb=1,wpro%nzb
       call weiR_write(X,izb)
    end do

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiR_assign
  !----------------------------------------------------------------
  subroutine weiR_report(X)
    type(weimag), 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)       ') "         *    Px: ",X%apx%n, X%apx%o, X%apx%d
    write(0,'(a,i3,f10.2,f10.3)       ') "         *    Py: ",X%apy%n, X%apy%o, X%apy%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)') "         *     Z: ",X%az%n,  X%az%o,  X%az%d," ->",wpro%nzb," 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 weiR_report
  !----------------------------------------------------------------
  subroutine weiR_read(X,k)
    type(weimag), intent(inout) :: X
    integer,      intent(in)    :: k
    integer :: stat
    character(len=128), parameter :: unit='weiR_read'
#ifdef DEBUG
    if(wpro%dbg>1) call in(unit)
#endif

    stat=sseek_block(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 weiR_read
  !----------------------------------------------------------------
  subroutine weiR_write(X,k)
    type(weimag), intent(inout) :: X
    integer,      intent(in)    :: k
    integer :: stat
    character(len=128), parameter :: unit='weiR_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 weiR_write
  !----------------------------------------------------------------
  subroutine weiR_mpicopy(X)
    type(weimag), intent(inout) :: X
    integer :: mpierr, msglen, msgsou
    integer :: izb
    integer :: nmx,nmy,nhx,nhy,nz
    character(len=128), parameter :: unit='weiR_mpicopy' 
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    nmx=size(X%bin,1)
    nmy=size(X%bin,2)
    nhx=size(X%bin,3)
    nhy=size(X%bin,4)
    nz =size(X%bin,5)

    msglen = nmx*nmy*nhx*nhy*nz

    msgsou=0
    do izb=1,wpro%nzb
       if(wpar%impi==0) write(0,'(a,a,i3)') weitab,"R send izb=",izb
       if(wpar%impi==0) call weiR_read(X,izb)
       call MPI_BARRIER(MPI_COMM_WORLD,mpierr)
       call MPI_BCAST(X%bin,msglen,MPI_COMPLEX,msgsou,MPI_COMM_WORLD,mpierr)
       if(wpar%impi/=0) call weiR_write(X,izb)
    end do

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weiR_mpicopy
  !----------------------------------------------------------------
  subroutine weiR_mpisumm(X)
    type(weimag), intent(inout) :: X
    integer :: mpierr, msglen, msgrec
    integer :: izb
    complex, pointer :: tmp(:,:,:,:,:)
    integer :: nmx,nmy,nhx,nhy,nz
    character(len=128), parameter :: unit='weiR_mpisumm' 
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif
    nmx=size(X%bin,1)
    nmy=size(X%bin,2)
    nhx=size(X%bin,3)
    nhy=size(X%bin,4)
    nz =size(X%bin,5)

    msglen = nmx*nmy*nhx*nhy*nz
    msgrec = 0

    if(wpar%impi==0) then
       call weimem_add(8,nmx*nmy*nhx*nhy*nz,unit)
       allocate(tmp(nmx,nmy,nhx,nhy,nz))
    else
       call weimem_add(8,1,unit)
       allocate(tmp(1,1,1,1,1))
    end if
    
    do izb=1,wpro%nzb
       if(wpar%impi==0) write(0,'(a,a,i3)') weitab,"R recv izb=",izb
       call weiR_read(X,izb)
       call MPI_BARRIER(MPI_COMM_WORLD,mpierr)
       call MPI_REDUCE(X%bin,tmp,msglen,MPI_COMPLEX,MPI_SUM,msgrec,MPI_COMM_WORLD,mpierr)
       if(wpar%impi==0) X%bin = tmp
       if(wpar%impi==0) call weiR_write(X,izb)
    end do

    if(wpar%impi==0) then
       call weimem_cut(8,nmx*nmy*nhx*nhy*nz,unit)
       deallocate(tmp)
    else
       call weimem_cut(8,1,unit)
       deallocate(tmp)
    end if

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

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

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






