module wei_obj_R
  use wei_util
  use wei_process
  use wei_type
  use wei_lagrange
  !---------------
  implicit none

  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!
  type(weimag), public :: R      ! background   image
  type(weimag), public ::dR      ! perturbation image
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!

contains
  !----------------------------------------------------------------
  subroutine weibimag_init(X,tag)
    type(weimag),       intent(inout) :: X
    character(len=128), intent(in)    :: tag

    X%tag = tag
    if(isfile(X%tag)) then
       call image_history_init(X)
    else
       call image_param_init(X)
    end if

    if(pro%use_disk_buf) then
       write(X%tag_var,"(a,a)") "R_var",pro%tag_num
       call auxinout(X%tag_var)
       call srite(X%tag_var,0.,4)
    end if

  end subroutine weibimag_init
  !----------------------------------------------------------------
  subroutine weidimag_init(X,tag)
    type(weimag),       intent(inout) :: X
    character(len=128), intent(in)    :: tag

    X%tag =tag
    if(isfile(X%tag)) then
       call image_history_init(X)
    else
       call image_param_init(X)
    end if

  end subroutine weidimag_init
  !----------------------------------------------------------------
  subroutine image_history_init(X)
    type(weimag), intent(inout) :: X
#ifdef DEBUG
    character(len=128), parameter :: unit='image_history_init'
    if(pro%debug) call in(unit)
#endif

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

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

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

       call hinaxis(X%tag,X%apy,4); X%apy%l='py'
       if(pro%apy%n ==0) then
          pro%apy = X%apy
       else
          if(.not. coax( X%apy, pro%apy)) then
             call image_report(X)
             call seperr('R.apy /= pro.apy')
          end if
       end if
    else
      call hinaxis(X%tag,X%ahx,3); X%ahx%l='hx'
      if(pro%ahx%n ==0) then
         pro%ahx = X%ahx
      else
         if(.not. coax( X%ahx, pro%ahx)) then
            call image_report(X)
            call seperr('R.ahx /= pro.ahx')
         end if
      end if

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

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


    call from_aux(trim(X%tag),"esize",X%esize,4)
    call image_setup(X)

#ifdef SEP_MPI
    if(pro%mpi_nthread>1) then
       X%tag_0=X%tag;
       write(X%tag,"(a,a)") trim(X%tag_0),pro%tag_num
       call auxscr(X%tag)
       call houaxis(X%tag,X%amx,1); call houaxis(X%tag,X%amy,2)
       call houaxis(X%tag,X%ahx,3); call houaxis(X%tag,X%ahy,4)
       call houaxis(X%tag,X%az, 5)
       call to_history("esize",X%esize,X%tag)
    end if
#endif

#ifdef DEBUG
    if(pro%debug) call image_report(X)
    if(pro%debug) call out(unit)
#endif
  end subroutine image_history_init
  !----------------------------------------------------------------
  subroutine image_param_init(X)
    type(weimag), intent(inout) :: X
    logical :: image_real
#ifdef DEBUG
    character(len=128), parameter :: unit='image_param_init'
    if(pro%debug) call in(unit)
#endif

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

    X%amx = pro%amx; call houaxis(X%tag,X%amx,1)
    X%amy = pro%amy; call houaxis(X%tag,X%amy,2)
    if(pro%pCIG) then
       X%apx = pro%apx; call houaxis(X%tag,X%apx,3)
       X%apy = pro%apy; call houaxis(X%tag,X%apy,4)
    else
       X%ahx = pro%ahx; call houaxis(X%tag,X%ahx,3)
       X%ahy = pro%ahy; call houaxis(X%tag,X%ahy,4)
    end if
    X%az  = pro%az;  call houaxis(X%tag,X%az, 5)

    call from_param("image_real",image_real,.false.)
    if(image_real) then
       X%esize=4
    else
       X%esize=8
    end if

    call houaxis(X%tag,X%amx,1);
    call houaxis(X%tag,X%amy,2)

    if(pro%pCIG) then
      call houaxis(X%tag,X%apx,3); call houaxis(X%tag,X%apy,4)
    else
		call houaxis(X%tag,X%ahx,3); call houaxis(X%tag,X%ahy,4)
    end if
    call houaxis(X%tag,X%az, 5)

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

    call image_setup(X)
    if(pro%mpi_ithread==0) then
       call image_assign(X,cmplx(0.))
    end if

    if(pro%mpi_nthread>1) then
       X%tag_0=X%tag;
       write(X%tag,"(a,a)") trim(X%tag_0),pro%tag_num
       call auxscr(X%tag)
       call houaxis(X%tag,X%amx,1);call houaxis(X%tag,X%amy,2)
       call houaxis(X%tag,X%ahx,3); call houaxis(X%tag,X%ahy,4)
       call houaxis(X%tag,X%az, 5)
       call to_history("esize",X%esize,X%tag)
    end if

#ifdef DEBUG
    if(pro%debug) call image_report(X)
    if(pro%debug) call out(unit)
#endif
  end subroutine image_param_init
  !--------------------------------------------------------------
  subroutine image_setup(X)
    type(weimag), intent(inout) :: X
    integer :: er
#ifdef DEBUG
    character(len=128), parameter :: unit='image_setup'
    if(pro%debug) call in(unit)
#endif

    !! DEPTH blocking
    if( (1.0*X%az%n)/pro%nzs /= X%az%n/pro%nzs ) then
       call image_report(R)
       write(0,*) "COMPARE",X%az%n,pro%nzs,X%az%n
       call seperr('R.az.n is not a multiple of nzs')
    end if

    pro%nzb = pro%az%n / pro%nzs
    !!--------------------
    allocate(X%bin(pro%amx%n,pro%amy%n,pro%ahx%n,pro%ahy%n,pro%nzs),stat=er)
    if(er/=0) call seperr('cant allocate R.bin')
    X%bin=(0.0,0.0)
    !!--------------------
    X%bites=X%esize*product((/pro%amx%n,pro%amy%n,pro%ahx%n,pro%ahy%n,pro%nzs/))

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine image_setup
  !----------------------------------------------------------------
  subroutine image_assign(X,v)
    type(weimag), intent(inout) :: X
    complex, intent(in) :: v
    integer :: izb
#ifdef DEBUG
    character(len=128), parameter :: unit='image_assign'
    if(pro%debug) call in(unit)
#endif

    X%bin=v
    do izb=1,pro%nzb
       call image_write(X,izb)
    end do

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine image_assign
  !----------------------------------------------------------------
  subroutine image_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," ->",pro%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"
       write(0,'(a,i1)')                       "         * esize=",X%esize
    end if
    write(0,*)'        ************************************************'
  end subroutine image_report
  !----------------------------------------------------------------
  subroutine image_read(X,k,iw)
    type(weimag), intent(inout) :: X
    integer, optional :: iw
    integer, intent(in) :: k
    integer :: stat,i_w
    real, allocatable :: temp_buf(:,:,:,:,:)
#ifdef DEBUG
    character(len=128), parameter :: unit='image_read'
    if(pro%debug) call in(unit)
#endif

    i_w=0; if(present(iw)) i_w=iw

    if(i_w==0) then 

      stat=sseek_block(X%tag,k-1,X%bites,0)
      if(X%esize==4) then
         allocate(temp_buf(size(X%bin,1),size(X%bin,2),size(X%bin,3),&
         size(X%bin,4), size(X%bin,5)),stat=stat)
         if(stat/=0) call erexit("trouble allocating temp image file")
         call sreed(trim(X%tag), temp_buf, X%bites)
         X%bin=cmplx(temp_buf,0.)
         deallocate(temp_buf)
      else
         call sreed(trim(X%tag), X%bin, X%bites)
      end if
    else
      call seperr("unimplented, varaible sampling image read")
    end if

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine image_read
  !----------------------------------------------------------------
  subroutine combine_image(X)
    type(weimag), intent(inout) :: X
    integer ::  ierr,iadd
    integer, external :: sseek
    
#ifdef SEP_MPI
    iadd=0;if(pro%add) iadd=1;
    if(pro%mpi_nthread>1 .and. pro%adj) then
       if(pro%mpi_ithread==0) ierr= sseek(trim(x%tag_0),0,0)
       ierr=sseek(trim(X%tag),0,0) 
       if(0/=mpi_sep_tag_sum(pro%mpi_ithread,X%tag,0,X%tag_0,X%bites/4,&
       pro%nzb,iadd)) call seperr("trouble combining image \n");
    end if
!    call aux_unlink(X%tag)
    if(pro%use_disk_buf) call aux_unlink(X%tag_var)
#endif
  end subroutine combine_image
  !----------------------------------------------------------------
  subroutine distribute_image(X)
    type(weimag), intent(inout) :: X

#ifdef SEP_MPI
    if(pro%mpi_nthread>1 .and. .not. pro%adj) then
        call sseek(trim(X%tag),0,0)
       if(pro%mpi_ithread==0) then
          call sseek(trim(X%tag_0),0,0)
        end if
       if(0/=mpi_sep_tag_bcast(0,X%tag_0,pro%mpi_ithread,X%tag,X%bites/4,&
       pro%nzb)) call seperr("trouble distributing image \n")
    end if
#endif
  end subroutine distribute_image
  !----------------------------------------------------------------
  subroutine image_write(X,k,sum,iw,nz)
    logical, optional :: sum
    logical           :: sum_it
    integer, optional :: iw,nz
    integer           :: i_w,n_z
    type(weimag), intent(inout) :: X
    integer, intent(in) :: k
    integer :: iz,i,izlast,izfirst,ib,stat
    real    :: z,zcur
    logical :: write_out,stored
    complex, allocatable :: tempc_buf(:,:,:,:)
    complex, allocatable :: temp2_buf(:,:,:,:)
#ifdef DEBUG
    character(len=128), parameter :: unit='image_write'
    if(pro%debug) call in(unit)
#endif
    
    sum_it=.false.
    if(present(sum)) sum_it=sum
    i_w=0; if(present(iw)) i_w=iw
    n_z=0; if(present(nz)) n_z=nz


    !Just doing a normal write
    if(i_w==0) then  
      iz=(k-1)*size(X%bin,5)
      do i=1,size(X%bin,5)
        iz=iz+1
        call image_write_iz(X,X%bin(:,:,:,:,i),iz,sum_it)
      end do
    else
      allocate(tempc_buf(size(X%bin,1),size(X%bin,2),size(X%bin,3),&
        size(X%bin,4))); 
      allocate(temp2_buf(size(X%bin,1),size(X%bin,2),size(X%bin,3),&
        size(X%bin,4))); 

      !
      !
      !First check to see  what we can write
      !
      !

      izlast=min(size(pro%z_samp(i_w)%dz),size(X%bin,5)*k) !last iz step done
      zcur=pro%z_samp(i_w)%zloc(izlast) ! depth we have gone done to




      z=pro%az%o;iz=1;
      do while(z<=zcur .and. iz <= pro%az%n) 
        if(.not. pro%filled(iz)) then
          write_out=.true.
          do i=1,size(pro%z_samp(i_w)%iloc,1) 
            if(pro%z_samp(i_w)%iloc(i,iz) > izlast .and.  &
               abs(pro%z_samp(i_w)%frac(i,iz)) > .001) &
                write_out=.false.
          end do 
          if(write_out) then
            tempc_buf=0.
            call lagrange_interp(.true.,i_w,k,iz,X%bin,tempc_buf,temp2_buf,X%tag_var)
            call image_write_iz(X,tempc_buf,iz,sum_it)
            pro%filled(iz)=.true.
          end if 
         end if
         z=z+pro%az%d;iz=iz+1
       end do

      !
      !
      ! Check to see if we can clear any incore or disk memory
      !
      !

      ! Find the smallest iz still needed
      izfirst=izlast
      do iz=1, size(pro%filled) 
        if(.not. pro%filled(iz)) then
          do i=1,size(pro%z_samp(i_w)%iloc,1)
            if(abs(pro%z_samp(i_w)%frac(i,iz))>.001) &
              izfirst=min(izfirst,pro%z_samp(i_w)%iloc(i,iz)) 
          end do
        end if 
      end do 

      do i=1,pro%ndisk_buf
        if(pro%disk_buf(i) < izfirst ) pro%disk_buf(i)=0
      end do
    
      do i=1,pro%nmem_buf
        if(pro%mem_buf(i) < izfirst ) pro%mem_buf(i)=0
      end do
    
      ! 
      ! 
      !store any newly created image components that are still needed
      ! 
      ! 

      iz=(k-1)*size(X%bin,5)
      do i=1,n_z
        iz=iz+1
        if(iz>=izfirst) then !we need to store this image compoenent
          stored=.false.;ib=1;
          do while(ib <= pro%nmem_buf .and. .not. stored) 
            if(pro%mem_buf(ib)==0) then
              pro%v_i_buf(:,:,:,:,ib)=X%bin(:,:,:,:,i)
              stored=.true.; pro%mem_buf(ib)=iz
              pro%z_samp(iw)%status(iz)=ib
            end if 
            ib=ib+1; 
          end do
          ib=1;
          do while(ib <= pro%ndisk_buf .and. .not. stored) 
            if(pro%disk_buf(ib)==0) then
              stat=sseek_block(X%tag_var,ib-1,size(tempc_buf)*8,0)
              call srite(X%tag_var, X%bin(:,:,:,:,i), size(tempc_buf)*8)
              stored=.true.; pro%disk_buf(ib)=iz
              pro%z_samp(iw)%status(iz)=-ib
            end if 
            ib=ib+1;
          end do
        end if 
      end do
      deallocate(tempc_buf,temp2_buf)
    end if
  end subroutine image_write

  subroutine image_write_iz(X,image,k,sum_it)
    logical           :: sum_it
    type(weimag), intent(inout) :: X
    integer, intent(in) :: k
    integer :: stat
    real, allocatable :: tempr_buf(:,:,:,:)
    complex, allocatable :: tempc_buf(:,:,:,:)
    complex :: image(:,:,:,:)
    
#ifdef DEBUG
    character(len=128), parameter :: unit='image_write_iz'
    if(pro%debug) call in(unit)
#endif
    
    if(sum_it) then
       stat=sseek_block(X%tag,k-1,X%bites/size(X%bin,5),0)
       if(X%esize==4) then
          allocate(tempr_buf(size(X%bin,1),size(X%bin,2),size(X%bin,3),& 
          size(X%bin,4)))
          call sreed(X%tag, tempr_buf,size(tempr_buf)*4)
          image=image+cmplx(tempr_buf,0.)
          deallocate(tempr_buf)
       else
          allocate(tempc_buf(size(X%bin,1),size(X%bin,2),size(X%bin,3),& 
          size(X%bin,4)))
          call sreed(X%tag, tempc_buf,size(tempc_buf)*8)
          image=image+tempc_buf
          deallocate(tempc_buf)
       end if
    end if
    
    stat=sseek_block(X%tag,k-1,X%bites/size(X%bin,5),0)
    if(X%esize==8) then
       call srite(X%tag, image, size(image)*8)
    else 
       call srite(X%tag, real(image), size(image)*4)
    end if

    call to_history("zblock",k,X%tag)
    
#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine image_write_iz
  !----------------------------------------------------------------
end module wei_obj_R

