program TransfSP
  use sep
  use flags
  use timers  
  implicit none

  integer :: i1,i2,i3,i4,i5,n1,n2,n3,n4,n5
  complex :: v
  integer*8,pointer:: planI(:)
  complex,  pointer:: cdata(:,:),cdata2(:,:)
  real,     pointer:: data(:,:)
  real::startFFT,finishFFT,startW,finishW,d1
  call sep_init(SOURCE)
  call from_history("n1",n1)
  call from_history("d1",d1)
  call from_history("n2",n2)
  call from_history("n3",n3)
  call from_history("n4",n4)
  call from_history("n5",n5)

  call from_param("f_min",f_min,0.)
  call from_param("f_max",f_max,0.5/dt)
  call from_param("f_min1",f_min1,f_min+df*2)
  call from_param("f_max1",f_max1,f_max-df*3)
  of=0.
  i_f_min=max(1,nint(f_min/df))
  i_f_max=min(int(nt/2),nint(f_max/df))
  nf=i_f_max-i_f_min+1
  f_min = (i_f_min-1)*df
  f_max = (i_f_max-1)*df

  call to_history("n1",n2);call to_history("label1",'x')
  call to_history("n2",n3);call to_history("label2",'y')
  call to_history("n3",nf);call to_history("label3",'f')
  call to_history("n4",n4);call to_history("label4",'sx')
  call to_history("n5",n5);call to_history("label5",'sy')

  call to_history("d1",d2);  call to_history("o1",o1)
  call to_history("d2",d3);  call to_history("o1",o2)
  call to_history("d3",df);  call to_history("o1",of)
  call to_history("d4",d4);  call to_history("o1",o4)
  call to_history("d5",d5);  call to_history("o1",o5)
  call to_history("esize",8)

  allocate(data(n1,n2),cdata(n1,n2),cdata2(n2,nf),planI(n2))
  write(0,*) 'planing'
  do i2=1,n2
     call sfftw_plan_dft_1d(planI(i2),n1,&
     cdata(1,i2),cdata(1,i2),FFTW_BACKWARD,FFTW_PATIENT)        
  end do
  sy:do i5=1,n5
     sx:do i4=1,n4
        gy:do i3=1,n3

           call sep_read(data)
!!$OMP PARALLEL DO PRIVATE( i2,i1,v)
           gx:do i2=1,n2
              t: do i1=1,n1
                 cdata(i1,i2)=cmplx(data(i1,i2),0.)
              end do t
              call sfftw_execute(planI(i2))
           end do gx
           tranx:do i2=1,n2
              tranw:do i1=1,nf
                 cdata2(i2,i1)=cdata(i1,i2)
              end do tranw
           end do tranx

              w:do i1=1,n1
                 v=cdata(i1,i2)
                 cdata2(i1,i2)=v
                 cdata2(i1+n1,i2)=v
              end do w
              call start_timer(startFFT)

              call end_timer(finishFFT,startFFT)
              call start_timer(startFFT)
              w2:do i1=1,n1
                 data(i1,i2)=real(cdata(i1,i2))
              end do w2
           end do gx
           call start_timer(startW)
           call sep_write(data)
           call end_timer(finishW,startW)
        end do gy
!!$OMP END PARALLEL DO 
     end do sx
  end do sy

write(0,'(8x,a,f10.3,a)')' *       Time FFTing       :', &
  finishFFT,' (s)'
write(0,'(8x,a,f10.3,a)')' *       Time Writing       :', &
  finishW,' (s)'
deallocate(planI,cdata,cdata2)
write(0,*) 'Done'
end program FToneside
