module  wei_file_mod
  use wei_wavespace_type_mod
  use sep3d_struct_mod
  use sep
  implicit none
  real, private, save :: mx

  type wfile
    integer            :: n(5)
    real               :: o(5) 
    real               :: d(5)
    character(len=128) :: label(5)
    character(len=128) :: tag
    logical            :: input
    logical            :: add
    logical            :: float
    type(wspace),pointer       :: wave
  end type


  contains







!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!     INIT ROUTINIES          !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine init_file_tag(fdes,wave,tag,float,freq,axes)
    type(wfile)           :: fdes
    type(wspace),pointer          :: wave
    character(len=*)      :: tag
    logical               :: float   !whether or not dataset is complex
    logical               :: freq    !whether or not this is a frequency file
    character(len=128)    :: axes(5) !names   of axes for error checking
    type(sep3d)           :: struct
    integer               :: i,n
    real                  :: o,d

#ifdef DEBUG 
    write(0,*) "IN INIT FILE TAG"
#endif

    call from_param("mx",mx,65536.)
    call init_sep3d(tag,struct,"INPUT")
!    if(sep3d_ndims(struct)/=5) call seperr("expecting 5-D dataset")

    do i=1,5
      if(i < 5)  then
        if(abs(struct%o(i)-wave%o(i)) > abs(struct%d(i))/5.) then
          write(0,*) "AXIS=",trim(axes(i)), " tag=",trim(tag)
          write(0,*) "TAG-0=",struct%o(i),"DES-o=",wave%o(i)
          call seperr("Os don't match")
        end if
        if(abs(struct%d(i)/wave%d(i)-1.) > abs(struct%d(i))/5.) then
          write(0,*) "AXIS=",trim(axes(i)), " tag=",trim(tag)
          write(0,*) "TAG-d=",struct%d(i),"DES-d=",wave%d(i)
          call seperr("ds don't match")
        end if
        if(struct%n(i)/=wave%n(i)) then
          write(0,*) "AXIS=",trim(axes(i)), " tag=",trim(tag)
          write(0,*) "TAG-n=",struct%n(i),"DES-n=",wave%n(i)
          call seperr("ns don't match")
        end if
      else
        if(  freq) then
          fdes%n(5)=size(wave%ws)
          fdes%o(5)=wave%ws(1)
          fdes%d(5)=wave%ws(2)-wave%ws(1)
          if(struct%n(i)/=size(wave%ws)) then
            write(0,*) "AXIS=",trim(axes(i)), " tag=",trim(tag)
            write(0,*) "TAG-n=",struct%n(i),"DES-n=",size(wave%ws)
            call seperr("ns don't match")
          end if
        end if
        if(.not. freq) then
          fdes%n(5)=size(wave%dz)
          if(struct%n(i)/=size(wave%dz)) then
            write(0,*) "AXIS=",trim(axes(i)), " tag=",trim(tag)
            write(0,*) "TAG-n=",struct%n(i),"DES-n=",size(wave%dz)
            call seperr("ns don't match")
          end if
        end if
      end if
    end do

    fdes%float =float
    if(struct%data_format(1:4)=="COMP" .and.  fdes%float) then
      call seperr("File is real, initialized requested complex ")
    end if

    if(struct%data_format(1:4)=="FLOA" .and. .not. fdes%float) then
      call seperr("File is complex, initialized requested real ")
    end if


    fdes%n(1:4)=wave%n(1:4)
    fdes%o(1:4)=wave%o(1:4)
    fdes%d(1:4)=wave%d(1:4)
    fdes%tag   =tag
    fdes%wave  =>wave
#ifdef DEBUG 
    write(0,*) "OUT INIT FILE TAG"
#endif
  end subroutine



  !INITIALIZE A FILE FROM WAVEFIELD INFORMATION
  subroutine init_file_struct(fdes,wave,tag,float,freq,axes)
    type(wfile)           :: fdes
    type(wspace),pointer          :: wave
    character(len=*)      :: tag
    logical               :: float   !whether or not dataset is complex
    logical               :: freq    !whether or not this is a frequency file
    character(len=128)    :: axes(5) !names   of axes for error checking
    type(sep3d)           :: struct
    integer               :: i

#ifdef DEBUG 
    write(0,*) "IN INIT FILE STRUCT"
#endif

    fdes%n(1:4)=wave%n(1:4)
    fdes%o(1:4)=wave%o(1:4)
    fdes%d(1:4)=wave%d(1:4)
    if(freq) then
      fdes%n(5)=size(wave%ws)
      fdes%o(5)=wave%ws(1)/8./atan(1.)
      fdes%d(5)=(wave%ws(2)-wave%ws(1))/8./atan(1.)
    else
      fdes%n(5)=size(wave%dz)
      fdes%o(5)=wave%z0
      fdes%d(5)=wave%dz(1)
    end if

    do i=1,5
      fdes%label(i)=axes(i)
    end do

     
    if(.not. float) then
      call init_sep3d(struct,"OUTPUT","COMPLEX","REGULAR",fdes%n,fdes%o,fdes%d,fdes%label)
    else
      call init_sep3d(struct,"OUTPUT","FLOAT","REGULAR",fdes%n,fdes%o,fdes%d,fdes%label)
    end if

    fdes%tag   =tag
    fdes%float =float
    fdes%wave =>wave
    call sep3d_write_description(fdes%tag,struct)

    call auxclose(fdes%tag)
    call auxinout(fdes%tag)

#ifdef DEBUG 
    write(0,*) "OUT INIT FILE STRUCT"
#endif
  end subroutine


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!     READ ROUTINIES          !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  !READ SLICE FROM DATA
  subroutine read_i5(fdes,ifile,ibuf,wfld)
    type(wfile)         :: fdes
    integer             :: ifile !i5 of file
    integer             :: ibuf  !i5 of buffer
    complex             :: wfld(:,:,:,:,:)  !buffer to read into
    integer             :: i1,i2,i3,i4

#ifdef DEBUG 
    write(0,*) "IN READ i5"
#endif
  
    do i4=1,size(wfld,4);do i3=1,size(wfld,3);do i2=1,size(wfld,2);do i1=1,size(wfld,1)
      wfld(i1,i2,i3,i4,ibuf)=0.
    end do; end do; end do;end do
    if(.not. fdes%float) then  
      call read_i5_c(fdes,ifile,ibuf,wfld)
    else
      call read_i5_r(fdes,ifile,ibuf,wfld)
    end if
   
#ifdef DEBUG 
    write(0,*) "OUT READ i5"
#endif
    
  end subroutine

  subroutine read_i5_r(fdes,ifile,ibuf,wfld)
    type(wfile)         :: fdes
    integer             :: ifile !i5 of file
    integer             :: ibuf  !i5 of buffer
    complex             :: wfld(:,:,:,:,:)  !buffer to read into
    real,allocatable    :: buf(:) !temporary buffer
    integer             :: i3,i4,i1,i2,i,ibcx,ibcy,ibhx,ibhy

#ifdef DEBUG 
    write(0,*) "IN READ i5 R"
#endif

    ibcx=(fdes%wave%npad(1)-fdes%wave%n(1))/2
    ibcy=(fdes%wave%npad(2)-fdes%wave%n(2))/2
    ibhx=(fdes%wave%npad(3)-fdes%wave%n(3))/2
    ibhy=(fdes%wave%npad(4)-fdes%wave%n(4))/2


    allocate(buf(fdes%n(1)*fdes%n(2)))
    do i4=1,fdes%n(4)
      do i3=1,fdes%n(3)
        call read_slice_r(fdes,ifile,i3,i4,ifile,buf)
        i=0
        do i2=1,fdes%n(2)
          do i1=1,fdes%n(1)
            i=i+1
            wfld(i1+ibcx,i2+ibcy,i3+ibhx,i4+ibhy,ibuf)=cmplx(buf(i),0.)
          end do
        end do
      end do
    end do
    deallocate(buf)
#ifdef DEBUG 
    write(0,*) "OUT READ i5 R"
#endif
  end subroutine

   
  subroutine read_i5_c(fdes,ifile,ibuf,wfld)
    type(wfile)         :: fdes
    integer             :: ifile !i5 of file
    integer             :: ibuf  !i5 of buffer
    complex             :: wfld(:,:,:,:,:)  !buffer to read into
    complex,allocatable :: buf(:) !temporary buffer
    integer             :: i3,i4,i1,i2,i,ibhx,ibhy,ibcx,ibcy

#ifdef DEBUG 
    write(0,*) "IN READ i5 C"
#endif

    ibcx=(fdes%wave%npad(1)-fdes%wave%n(1))/2
    ibcy=(fdes%wave%npad(2)-fdes%wave%n(2))/2
    ibhx=(fdes%wave%npad(3)-fdes%wave%n(3))/2
    ibhy=(fdes%wave%npad(4)-fdes%wave%n(4))/2


    allocate(buf(fdes%n(1)*fdes%n(2)))
    do i4=1,fdes%n(4)
      do i3=1,fdes%n(3)
        call read_slice_c(fdes,ifile,i3,i4,ifile,buf)
        i=0
        do i2=1,fdes%n(2)
          do i1=1,fdes%n(1)
            i=i+1
            wfld(i1+ibcx,i2+ibcy,i3+ibhx,i4+ibhy,ibuf)=buf(i)
          end do
        end do
      end do
    end do
    deallocate(buf)
#ifdef DEBUG 
    write(0,*) "OUT READ i5 C"
#endif
  end subroutine

  subroutine read_slice_r(fdes,ifile,i3,i4,i5,buf)
    type(wfile)       :: fdes
    integer           :: ifile
    integer           :: i3,i4,i5
    real              :: buf(:)
    integer,external  :: sreed_window

    
#ifdef DEBUG 
    write(0,*) "IN READ SLICE R"
#endif
    if(0/=sreed_window(fdes%tag,5,fdes%n,(/fdes%n(1),fdes%n(2),1,1,1/),(/0,0,i3-1,i4-1,i5-1/), &
      (/1,1,1,1,1/),4,buf))  then
       write(0,*) "trouble reading from ",trim(fdes%tag)
       call seperr("")
    end if
#ifdef DEBUG 
    write(0,*) "OUT READ SLICE R"
#endif

  end subroutine


  subroutine read_slice_c(fdes,ifile,i3,i4,i5,buf)
    type(wfile)       :: fdes
    integer           :: ifile
    integer           :: i3,i4,i5
    complex           :: buf(:)
    integer,external  :: sreed_window

    
#ifdef DEBUG 
    write(0,*) "IN READ SLICE C"
#endif
    if(0/=sreed_window(fdes%tag,5,fdes%n,(/fdes%n(1),fdes%n(2),1,1,1/),(/0,0,i3-1,i4-1,i5-1/), &
      (/1,1,1,1,1/),8,buf))  then
       write(0,*) "trouble reading from ",trim(fdes%tag)
       call seperr("")
    end if
#ifdef DEBUG 
    write(0,*) "OUT READ SLICE C"
#endif

  end subroutine


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!     WRITE ROUTINIES          !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


  subroutine write_i5(fdes,ifile,ibuf,wfld)
    type(wfile)         :: fdes
    integer             :: ifile !i5 of file
    integer             :: ibuf  !i5 of buffer
    complex             :: wfld(:,:,:,:,:)  !buffer to read into

#ifdef DEBUG 
    write(0,*) "IN WRITE i5"
#endif
    
    if(.not. fdes%float) then 
      call write_i5_c(fdes,ifile,ibuf,wfld)
    else
      call write_i5_r(fdes,ifile,ibuf,wfld)
    end if

#ifdef DEBUG 
    write(0,*) "OUT WRITE i5"
#endif
  end subroutine



  subroutine write_i5_c(fdes,ifile,ibuf,wfld)
    type(wfile)         :: fdes
    integer             :: ifile !i5 of file
    integer             :: ibuf  !i5 of buffer
    complex             :: wfld(:,:,:,:,:)  !buffer to read into
    complex,allocatable :: buf(:) !temporary buffer
    complex,allocatable :: buf2(:) !temporary buffer
    integer             :: i3,i4,i1,i2,i,ibcx,ibcy,ibhx,ibhy

#ifdef DEBUG 
    write(0,*) "IN WRITE i5 C"
#endif


    ibcx=(fdes%wave%npad(1)-fdes%wave%n(1))/2
    ibcy=(fdes%wave%npad(2)-fdes%wave%n(2))/2
    ibhx=(fdes%wave%npad(3)-fdes%wave%n(3))/2
    ibhy=(fdes%wave%npad(4)-fdes%wave%n(4))/2

    allocate(buf(fdes%n(1)*fdes%n(2)))
    if(fdes%add) allocate(buf2(fdes%n(1)*fdes%n(2)))


    do i4=1,fdes%n(4)
      do i3=1,fdes%n(3)
        if(fdes%add) then
          call read_slice_c(fdes,ifile,i3,i4,ifile,buf2)
          i=0
          do i2=1,fdes%n(2)
            do i1=1,fdes%n(1)
              i=i+1
              buf(i)=(wfld(i1+ibcx,i2+ibcy,i3+ibhx,i4+ibhy,ibuf))+buf2(i)
            end do
          end do
        else
          i=0
          do i2=1,fdes%n(2)
            do i1=1,fdes%n(1)
              i=i+1
              buf(i)=(wfld(i1+ibcx,i2+ibcy,i3+ibhx,i4+ibhy,ibuf))
            end do
          end do
        end if
        call write_slice_c(fdes,ifile,i3,i4,ifile,buf)
      end do
    end do

   deallocate(buf)
   if(fdes%add) deallocate(buf2)

#ifdef DEBUG 
    write(0,*) "OUT WRITE i5 C"
#endif

  end subroutine

  subroutine write_i5_r(fdes,ifile,ibuf,wfld)
    type(wfile)         :: fdes
    integer             :: ifile !i5 of file
    integer             :: ibuf  !i5 of buffer
    complex             :: wfld(:,:,:,:,:)  !buffer to read into
    real,allocatable    :: buf(:) !temporary buffer
    real,allocatable    :: buf2(:) !temporary buffer
    integer             :: i3,i4,i1,i2,i,ibcx,ibcy,ibhx,ibhy

#ifdef DEBUG 
    write(0,*) "IN WRITE i5 R"
#endif

    ibcx=(fdes%wave%npad(1)-fdes%wave%n(1))/2
    ibcy=(fdes%wave%npad(2)-fdes%wave%n(2))/2
    ibhx=(fdes%wave%npad(3)-fdes%wave%n(3))/2
    ibhy=(fdes%wave%npad(4)-fdes%wave%n(4))/2


    allocate(buf(fdes%n(1)*fdes%n(2)))
    if(fdes%add) allocate(buf2(fdes%n(1)*fdes%n(2)))
!    if(ifile > 179) call srite("itest",(wfld(181:600,:,:,:,ibuf)),size(wfld)/size(wfld,5)/size(wfld,1)*420*8)
    call srite("itest",wfld(:,:,:,:,ibuf),size(wfld)/size(wfld,5)*8);
;
;
    do i4=1,fdes%n(4);
      do i3=1,fdes%n(3);
        if(fdes%add) then;
          call read_slice_r(fdes,ifile,i3,i4,ifile,buf2)
          i=0
          do i2=1,fdes%n(2)
            do i1=1,fdes%n(1)
              i=i+1
              buf(i)=real(wfld(i1+ibcx,i2+ibcy,i3+ibhx,i4+ibhy,ibuf))++buf2(i)
            end do
          end do
!            if(i3==45) write(0,*) "IN",sum(buf2),sum(wfld(:,:,45,:,:))
        else
          i=0
          do i2=1,fdes%n(2)
            do i1=1,fdes%n(1)
              i=i+1

              buf(i)=real(wfld(i1+ibcx,i2+ibcy,i3+ibhx,i4+ibhy,ibuf))
            end do
          end do
        end if
        call write_slice_r(fdes,ifile,i3,i4,ifile,buf)
      end do
    end do

   deallocate(buf)
   if(fdes%add) deallocate(buf2)
#ifdef DEBUG 
    write(0,*) "OUT WRITE i5 R"
#endif

  end subroutine

  subroutine write_slice_c(fdes,ifile,i3,i4,i5,buf)
    type(wfile)       :: fdes
    integer           :: ifile
    integer           :: i3,i4,i5
    complex           :: buf(:)
    integer,external  :: srite_window
#ifdef DEBUG 
    write(0,*) "IN WRITE SLICE C"
#endif

    
    if(0/=srite_window(fdes%tag,5,fdes%n,(/fdes%n(1),fdes%n(2),1,1,1/),(/0,0,i3-1,i4-1,i5-1/), &
      (/1,1,1,1,1/),8,buf))  then
       write(0,*) "trouble reading from ",trim(fdes%tag)
       call seperr("")
    end if
#ifdef DEBUG 
    write(0,*) "OUT WRITE SLICE C"
#endif

  end subroutine





  subroutine write_slice_r(fdes,ifile,i3,i4,i5,buf)
    type(wfile)       :: fdes
    integer           :: ifile
    integer           :: i3,i4,i5
    real              :: buf(:)
    integer,external  :: srite_window
#ifdef DEBUG 
    write(0,*) "IN WRITE SLICE R"
#endif

    
    if(0/=srite_window(fdes%tag,5,fdes%n,(/fdes%n(1),fdes%n(2),1,1,1/),(/0,0,i3-1,i4-1,i5-1/), &
      (/1,1,1,1,1/),4,buf))  then
       write(0,*) "trouble reading from ",trim(fdes%tag)
       call seperr("")
    end if
#ifdef DEBUG 
    write(0,*) "OUT WRITE SLICE R"
#endif

  end subroutine






   


end module