module R_mod

  use util_mod
  use external_mod
  use down_mod

  implicit none
  integer, private, save,allocatable  :: image_master(:)
  integer, private, save  :: image_comm_group
  logical, private, save  :: image_lead_thread
  real,allocatable,dimension(:,:)::tmpr
  complex,allocatable,dimension(:,:):: tmpc
  logical,private,save::notall=.true.

contains
!----------------------------------------------------------------
! Image initializing routine
  logical function  image_init(X,tag)
    type(imag),       intent(inout) :: X
    character(len=128), intent(in)    :: tag
    integer ::ierr
    real::a
#ifdef DEBUG
    character(len=128), parameter :: unit='image_init'
    if(run%debug) call in(unit)
#endif
    image_init=F
    X%tag = tag

    if( down%add .and. isfile(X%tag) ) then
       if(.not. image_history_init(X)) then
          write(run%io,*) "trouble initializing image_history"
          return
       end if
    else
       if(.not. image_param_init(X)) then  !THIS ONE
          write(run%io,*) "trouble initializing image_param"
          return
       end if
    end if

    if(down%image_real) then
       allocate(tmpr(down%amx%n,down%amy%n),stat=ierr)
       if(ierr/=0) call seperr('cant allocate tmpr R.F90',ierr)
       tmpr=0.
    else
       allocate(tmpc(down%amx%n,down%amy%n),stat=ierr)
       if(ierr/=0) call seperr('cant allocate tmpc R.F90',ierr)
       tmpc=cmplx(0.,0.)
    end if

    a=4*size(tmpr)+4*size(tmpc)
    down%mem=down%mem+a
#ifdef MEMCHECK
    if(run%debug) write(run%io,*) 'R.F90 tmpr tmpc: allocate',a/1024./1024.,'Mbytes.  Total=',down%mem/1024./1024.
#endif
    image_init=T
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end function image_init

!----------------------------------------------------------------
! Image initializing routine
  logical function  aux_image_init(X,Y,tag)
    type(imag)       :: X
    type(imag),       intent(in) :: Y
    character(len=128), intent(in)    :: tag
    integer ::ierr
#ifdef DEBUG
    character(len=128), parameter :: unit='image_init'
    if(run%debug) call in(unit)
#endif
    aux_image_init=F
    X%tag=tag

    if( down%add .and. isfile(X%tag) ) then
       if(.not. image_history_init(X)) then
         write(run%io,*) "trouble initializing aux_image_history"
         return
       end if
       if(X%amx%n.ne.Y%amx%n.or.X%amy%n.ne.Y%amy%n.or.X%az%n.ne.Y%az%n) &
         call seperr('aux image file not the same size as R')
    else !!for new files    
       X%esize=Y%esize ;X%amx=Y%amx; X%ahx=Y%ahx; X%amy=Y%amy; X%ahy=Y%ahy; X%az=Y%az; X%aw=Y%aw
       X%indx=Y%indx; X%indy=Y%indy; X%finx=Y%finx; X%finy=Y%finy; X%from=Y%from; X%bites=Y%bites
       if(.not. create_image_file(X,X%tag)) then
          write(run%io,*) "trouble creating image"
          return
       end if
       if(.not. image_setup(X,.not.down%add)) then
          write(run%io,*) "trouble in image setup"
          return
       end if
       call image_report(X)
    end if

    aux_image_init=T
#ifdef DEBUG
      if(run%debug) call out(unit)
#endif
  end function aux_image_init

!----------------------------------------------------------------
! Initialize an image file based on parameters
  logical function image_param_init(X)
    type(imag), intent(inout) :: X
    integer :: nhx,nhx2,nhy,nhy2,noffx,noffy
    real::o
    logical:: m
#ifdef DEBUG
    character(len=128), parameter :: unit='image_param_init'
    if(run%debug) call in(unit)
#endif
    image_param_init=F
    X%from = "parameter"
    X%esize=4;if(.not.down%image_real) X%esize=8
!! . . depth o-rama
    X%az%d=down%az%d
    X%az%n=down%az%n-down%firstz+1
    X%az%o=down%az%o + (down%firstz-1)*down%az%d

!! . . calculate everything for CIGs
    call from_param('ahx_n',nhx,1)!! old style
    call from_param('ahy_n',nhy,1)
    call from_param('nhx',nhx2,1) !! new style
    call from_param('nhy',nhy2,1)
    down%ahx%n=max(nhx,nhx2)
    down%ahy%n=max(nhy,nhy2)

    call from_param('jhx',down%ahx%j,1)
    call from_param('jhy',down%ahy%j,1)
    call from_param('dhx',down%ahx%d,down%amx%d*down%ahx%j)
    call from_param('dhy',down%ahy%d,down%amy%d*down%ahy%j)

    if(down%ahx%n.gt.1) then
       !! make sure dhx/dx = integer
       nhx=nint(down%ahx%d/down%amx%d)
write(0,*)'dhx=',down%ahx%d
write(0,*)'dx=',down%amx%d
       o=nhx*down%amx%d
write(0,*)'nhx=',nhx
write(0,*)'o=',o
       !if(o.ne.down%ahx%d) then
       if(abs(o-down%ahx%d)>0.01) then
          write(run%io,"(a,f10.4)") " ** Supplied dhx not multiple of dx. Changed to",o
          down%ahx%d=o
       end if
       !! should have odd number of offsets as zero is included
       if(mod(down%ahx%n,2)==0) down%ahx%n=down%ahx%n+1 

       if(mod(int(10000.*down%ahx%d),int(10000.*down%amx%d)).ne.0) call seperr("dhx must be multiple of dx")
       o= -(real( down%ahx%n -1)/2.) * down%ahx%d
       call from_param("ohx",down%ahx%o,o)
    end if

    if(down%ahy%n.gt.1) then
       !! make sure dhx/dx = integer
       nhy=nint(down%ahy%d/down%amy%d)
write(0,*)'dhy=',down%ahy%d
write(0,*)'dy=',down%amy%d
       o=nhy*down%amy%d
write(0,*)'nhy=',nhy
write(0,*)'oy=',o
       !if(o.ne.down%ahy%d) then
       if(abs(o-down%ahy%d)>0.01) then
          write(run%io,"(a,f10.4)") " ** Supplied dhy not multiple of dy. Changed to",o
          down%ahy%d=o
       end if
       !! should have odd number of offsets as zero is included
       if(mod(down%ahy%n,2)==0) down%ahy%n=down%ahy%n+1 
       o= -(real( down%ahy%n -1)/2.) * down%ahy%d !offset origin begins w/ neg. offset
       call from_param("ohy",down%ahy%o,o)
    end if

    call from_param("hpos",m,.false.)
    if(m) then
       down%ahx%o=0.;down%ahy%o=0.
    end if
    call from_param("hneg",m,.false.)
    if(m) then
       down%ahx%o=-( down%ahx%n -1)* down%ahx%d
       down%ahy%o=-( down%ahy%n -1)* down%ahy%d
    end if

    X%ahx=down%ahx
    X%ahy=down%ahy

!!total image size from earth covered by recievers
    X%amx=down%aix 
    X%amy%n=1;X%amy%o=0.;X%amy%d=1.
    if (down%amy%n.gt.1)X%amy=down%aiy

    if(down%image_freqs) X%aw=down%aw; X%az%l='z'

    if(.not.down%add) then 
      if(.not. create_image_file(X,X%tag)) then
        write(run%io,*) "trouble creating image"
        return
       end if
    end if
    if(.not. image_setup(X,.not.down%add)) then
      write(run%io,*) "trouble in image setup"
      return
    end if

    call image_report(X)
    image_param_init=T
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end function image_param_init
!----------------------------------------------------------------
! Write out parameters for an image file
  logical function create_image_file(X,tag)
    type(imag) :: X
    type(dat)::D,D_d
    type(slow)::S
    type(twfld)::W,W_d
    logical :: i_real
    integer :: i,j,k
    character(len=*) ::  tag 
    create_image_file=F

    call auxscr(tag)
    call to_history("esize",X%esize,tag)
    if (down%image_freqs)then
       call houaxes(tag,X%amx,X%amy,X%ahx,X%ahy,X%aw,X%az)
    else
       call houaxes(tag,X%amx,X%amy,X%ahx,X%ahy,X%az)
       call to_history("n6",1 ,tag);call to_history("o6",1 ,tag)
       call to_history("d6",1 ,tag) !!sometimes you want some ones
    end if

    X%indx=1;X%indy=1
    i=srite(tag,0.,4)
    i=sseek(tag,0,0)

    create_image_file=T
  end function create_image_file
!----------------------------------------------------------------
  subroutine param_write(tag)
    character(len=*) ::  tag 

    !! General info into
    call to_history("nref",down%nref ,tag)
    call to_history("rollalong",down%roll,tag)
    call to_history("offset",down%offset,tag)
    call to_history("passive",down%passive ,tag)
    call to_history("fwdscat",down%fwdscat,tag)
    call to_history("Greens functions", down%datum,tag)
    call to_history("watercolumn",down%awtr%o,tag)
    call to_history("wtr_steps",down%awtr%n,tag)
    call to_history("velocity",down%vel ,tag)
    call to_history("operation",run%operation ,tag)
    call to_history("adj",down%adj ,tag)
    call to_history("add", down%add,tag)
    call to_history("cwave",down%cwave ,tag)
    call to_history("image_freqs",down%image_freqs,tag)
    call to_history("firstz",down%firstz,tag)
    call to_history("firstzW",down%firstzW,tag)
    call to_history("antlialiasing",down%antion,tag)
    call to_history("zblocks",down%bz%n)
    call to_history("zblock",down%zblock)
    call to_history("wblocks",down%bw%n)
    call to_history("wperblock",down%nws)
    if (down%antion) then
       call to_history("pmax_x",down%pmax_x,tag)
       call to_history("pmax_y",down%pmax_y,tag)
       call to_history("pmin_x",down%pmin_x,tag)
       call to_history("pmin_y",down%pmin_y,tag)
    end if
    call to_history("eps",down%eps,tag)
    call to_history("maxang",down%maxang,tag)
    call to_history("nx_shots",down%nx_shots ,tag)
    call to_history("ox_shots",down%ox_shots ,tag)
    call to_history("ny_shots",down%ny_shots ,tag)
    call to_history("oy_shots",down%oy_shots ,tag)
    call to_history("taper",down%taperme ,tag)
    call to_history("minx1",down%minx1,tag)
    call to_history("maxx1",down%maxx1,tag)
    call to_history("miny1",down%miny1,tag)
    call to_history("maxy1",down%maxy1,tag)
    call to_history("decon",down%dcnv ,tag)
    call to_history("smooth",down%smoothme ,tag)
    call to_history("rectx",down%rect1 ,tag)
    call to_history("recty",down%rect2 ,tag)
    call to_history("WIOu",down%WIOu,tag)
    call to_history("WIOd",down%WIOd,tag)
    call to_history("image",down%image,tag)
    call to_history("mult",down%mult,tag)
    call to_history("normd",down%normd,tag)
    call to_history("gathers",down%gathers,tag)
    call to_history("sourcefn",down%sourcefn,tag)
    call to_history("headD",down%headD,tag)
    call to_history("headD_d",down%headD_d,tag)
    call to_history("xyreverse",down%xyr,tag)
    call to_history("badFFTsign",down%badsign,tag)
    call to_history("nsx",down%asx%n)
    call to_history("osx",down%asx%o)
    call to_history("dsx",down%asx%d)
    call to_history("nsy",down%asy%n)
    call to_history("osy",down%asy%o)
    call to_history("dsy",down%asy%d)
    call to_history("totalshots",down%evnt%n)
  end subroutine param_write
!--------------------------------------------------------------
!  Core routine to setup an image
  logical function image_setup(X,zero)
    type(imag), intent(inout) :: X
    integer :: er,i,j,l
    real::a
    logical :: zero
    complex,allocatable,dimension(:,:)::tmp2
    real,allocatable,dimension(:,:)::tmp1
    integer:: i1,i2,i3,i4,i5

#ifdef DEBUG
    character(len=128), parameter :: unit='image_setup'
    if(run%debug) call in(unit)
#endif
    image_setup=F
    !! write a bunch of zeros
    if(zero) then
       if(X%esize==4) then
          allocate(tmp1(down%aix%n,down%aiy%n))
          tmp1=0.
          do i=1,X%az%n
             do j=1,down%ahy%n
                do l=1,down%ahx%n
                   er=srite(X%tag,tmp1,size(tmp1)*4)
                end do
             end do
          end do
          deallocate(tmp1)
       else
          allocate(tmp2(down%aix%n,down%aiy%n))
          tmp2=cmplx(0.,0.)
          do i=1,X%az%n
             do j=1,down%ahy%n
                do l=1,down%ahx%n
                   er=srite(X%tag,tmp2,size(tmp2)*8)
                end do
             end do
          end do
          deallocate(tmp2)
       end if
       er=sseek(X%tag,0,0)
    end if

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

    er=0
    if (.not. down%image_freqs) then
       allocate(X%bin(down%amx%n,down%amy%n,down%ahx%n,down%ahy%n,down%nzs),stat=er)
       if(er/=0) call seperr('cant allocate image sized binary',X%tag)
       do i5=1,down%nzs
          do i4=1,down%ahy%n
             do i3=1,down%ahx%n
                do i2=1,down%ay%n
                   do i1=1,down%ax%n
                      X%bin(i1,i2,i3,i4,i5)=cmplx(0.,0.)
                   end do
                end do
             end do
          end do
       end do
    else
       X%bites=X%esize*&
       product((/down%amx%n,down%amy%n,down%ahx%n,down%ahy%n,down%nws/))
       allocate(X%bin(down%amx%n,down%amy%n,down%ahx%n,down%ahy%n,1),stat=er)
       if(er/=0) call seperr('cant allocate image sized binary',X%tag)
       do i4=1,down%ahy%n
          do i3=1,down%ahx%n
             do i2=1,down%ay%n
                do i1=1,down%ax%n
                   X%bin(i1,i2,i3,i4,1)=cmplx(0.,0.)
                end do
             end do
          end do
       end do
    end if

    a=8*size(X%bin)
    down%mem=down%mem+a
#ifdef MEMCHECK
       if(run%debug) write(run%io,*) 'R.F90 bin: allocate',a/1024./1024.,'Mbytes.  Total=',down%mem/1024./1024.
#endif
    image_setup=T
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end function image_setup
!--------------------------------------------------------------
!  Assign value to the image file
  subroutine image_assign(X,v)
    type(imag), intent(inout) :: X
    complex, intent(in) :: v
    integer :: izb,n,i,ierr
    real:: u
#ifdef DEBUG
    character(len=128), parameter :: unit='image_assign'
    if(run%debug) call in(unit)
#endif
    u=real(v)
    X%bin=v
    ierr=sseek(X%tag,0,0)
    if (X%esize==4) then
       tmpr=u
       do i=1,down%az%n
          ierr=srite(X%tag,tmpr,size(tmpr)*4)
       end do
    else
       ierr=srite(X%tag,X%bin,size(X%bin)*8)
    end if

    ierr=sseek(X%tag,0,0)

#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end subroutine image_assign
!----------------------------------------------------------------
! Write info about image
  subroutine image_report(X)
    type(imag), intent(in) :: X
    real::bb
    bb=real(down%az%n)/real(down%nzs)
    write(run%io,*)'        ************************************************'
    write(run%io,*)'        *   ',trim(X%tag),' from ',trim(X%from)
    write(run%io,*)'        *                 [n]      [o]      [d]'
    write(run%io,'(a,i4,f10.3,f10.3)       ') "         *             x: ",X%amx%n, X%amx%o, X%amx%d
    write(run%io,'(a,i4,f10.3,f10.3)       ') "         *             y: ",X%amy%n, X%amy%o, X%amy%d
    write(run%io,'(a,i4,f10.3,f10.3)       ') "         *     Wrk areax: ",down%amx%n, down%amx%o, down%amx%d
    write(run%io,'(a,i4,f10.3,f10.3)       ') "         *     Wrk areay: ",down%amy%n, down%amy%o, down%amy%d
    write(run%io,'(a,i4,f10.3,f10.3)       ') "         *            hx: ",X%ahx%n, X%ahx%o, X%ahx%d
    write(run%io,'(a,i4,f10.3,f10.3)       ') "         *            hy: ",X%ahy%n, X%ahy%o, X%ahy%d
    write(run%io,'(a,i4,f10.3,f10.3,a,f5.2,a)') "         *             z: ",X%az%n,  X%az%o,  X%az%d," ->",bb," blocks"
    if(down%image_freqs)&
    write(run%io,'(a,i4,f10.3,f10.3)       ') "         *         freqs: ",X%aw%n, X%aw%o, X%aw%d
    write(run%io,'(a,f6.2,a)        ') "         *         block: ",2*X%bites/1024./1024.," Mbytes"
    write(run%io,'(a,i1)')                    "         *      esize= ",X%esize
  end subroutine image_report
!----------------------------------------------------------------
! Places single-shot image subspace into 
!          total image space. By brad 11/02/03
  subroutine image_location(X)
    type(imag), intent(inout)   ::X
    integer:: ierr,i,j,k,l,m,n,p,o
#ifdef DEBUG
    character(len=128), parameter :: unit='image_location'
    if(run%debug) call in(unit)
#endif
    X%indx= (1.5+ (down%ax%o-down%aix%o)/down%ax%d)
    X%indx= min(X%indx,down%aix%n-down%ax%n) !! Dont go off min/max points
    X%indx= max(X%indx,1)

    X%indy= (1.5+ (down%ay%o-down%aiy%o)/down%ay%d)
    X%indy= min(X%indy,down%aiy%n-down%ay%n) !! Dont go off min/max points
    X%indy= max(X%indy,1)

    X%finx=X%indx + down%ax%n -1
    X%finy=X%indy + down%ay%n -1
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end subroutine image_location
!----------------------------------------------------------------
! Function to read image
  subroutine image_read(X,k,zz)
    type(imag), intent(inout) :: X
    integer, intent(in) :: k,zz
    integer:: i1,i2,i3,i4,i5,ierr,n_dim_cube,n,m
    integer,dimension(5)::n_cube,n_wind,f_wind,j_wind
#ifdef DEBUG
    character(len=128), parameter :: unit='image_read'
    if(run%debug) call in(unit)
#endif
    write(run%io,*) "        *       Read image space ",trim(X%tag)

    n=zz

    n_dim_cube=5
    n_cube=(/   X%amx%n,   X%amy%n,X%ahx%n,X%ahy%n,X%az%n /)
    j_wind=(/1         ,1         ,1      ,1      ,1      /)
    n_wind=(/down%ax%n,down%ay%n  ,1      ,1      ,1      /)

    if(X%esize==4) then
       do i5=1,n          
          do i4=1,X%ahy%n
             do i3=1,X%ahx%n
                f_wind=(/X%indx -1,X%indy -1,i3-1,i4-1,i5-1+k /)
                ierr=sreed_window(X%tag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,X%esize,tmpr(1,1))
                if (ierr.lt.0) then
                   write(run%io,*) 'sreed_window return value=',ierr
                   call seperr("Baddness in image_read (real)")
                end if
                do i2=1,down%ay%n
                   do i1=1,down%ax%n
                      X%bin(i1,i2,i3,i4,i5)=cmplx(tmpr(i1,i2),0.)
                   end do
                end do
             end do
          end do
       end do
    else
       do i5=1,n          
          do i4=1,X%ahy%n
             do i3=1,X%ahx%n
                f_wind=(/X%indx -1,X%indy -1,i3-1,i4-1,i5-1+k /)
                ierr=sreed_window(X%tag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,X%esize,tmpc(1,1))
                if (ierr.lt.0) then
                   write(run%io,*) 'sreed_window return value=',ierr
                   call seperr("Baddness in image_read (complex)")
                end if

                do i2=1,down%ay%n
                   do i1=1,down%ax%n
                      X%bin(i1,i2,i3,i4,i5)=tmpc(i1,i2)
                   end do
                end do
             end do
          end do
       end do
    end if

#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end subroutine image_read
!----------------------------------------------------------------
! Image initializing routine when the image file exists
  logical function image_history_init(X)
    type(imag), intent(inout) :: X

#ifdef DEBUG
    character(len=128), parameter :: unit='image_history_init'
    if(run%debug) call in(unit)
#endif
    image_history_init=F

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

    call from_aux(X%tag,"esize",X%esize,4)
    call hinaxis(X%tag,X%amx,1); X%amx%l='x_'
    call hinaxis(X%tag,X%amy,2); X%amy%l='y_'
    call hinaxis(X%tag,X%az ,5); X%az%l ='z_'
    call hinaxis(X%tag,X%ahx,3); X%ahx%l='hx'
    call hinaxis(X%tag,X%ahy,4); X%ahy%l='hy'
    down%aix=X%amx
    down%aiy=X%amy
    down%ahx=X%ahx
    down%ahy=X%ahy

    if(.not. image_setup(X,.false.)) then
       write(run%io,*) "trouble setting up image"
       return
    end if
    call image_report(X)
    image_history_init=T

#ifdef DEBUG
    if(run%debug) call image_report(X)
    if(run%debug) call out(unit)
#endif
  end function image_history_init
!----------------------------------------------------------------
  subroutine update_index(X,Y)
    type(imag),intent(in)::Y
    type(imag),intent(inout)::X
    X%indx=Y%indx; X%indy=Y%indy; X%finx=Y%finx; X%finy=Y%finy;
  end subroutine update_index
!----------------------------------------------------------------
  subroutine imagearrays()
    integer::ierr
    if(down%image_real)then
       deallocate(tmpr)
       allocate(tmpr(down%ax%n,down%ay%n),stat=ierr)
       if(ierr/=0) call seperr('cant reallocate tmpr R.F90',ierr)
    else 
       deallocate(tmpc)
       allocate(tmpc(down%ax%n,down%ay%n),stat=ierr)
       if(ierr/=0) call seperr('cant reallocate tmpr R.F90',ierr)
    end if

  end subroutine imagearrays
!----------------------------------------------------------------
end module R_mod

