module cigrot

 use fftwrot
 use parmrot
 use typerot

 implicit none

 contains

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! subroutine sogrot_init(npadz_in,npadx_in,npadh_in)
!
! integer			:: npadz_in,npadx_in,npadh_in
!
! npadz=npadz_in; npadx=npadx_in; npadh=npadh_in
!
! end subroutine sogrot_init 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine sogrot(incig,outcig)

 real				:: incig(:,:,:),outcig(:,:,:)
 complex,allocatable		:: Cin(:,:,:),Cout(:,:,:),tfld3(:,:,:),tfld2(:,:),tfld1(:),temp(:),krow(:)
 integer				:: ix, iz, ih, ikh, nodes,k1,k2
 real				:: ta,tg,tteta,khrot,f,g
 real,allocatable			:: shift(:)

! if (node.eq.9999) then
!   !$OMP PARALLEL
!   nodes = omp_get_num_threads()
!   !$OMP END PARALLEL
! else
!   nodes=node
! end if
! write(0,*) "****************"
! write(0,*) "* NUM.THREADS = ",nodes
! write(0,*) "****************"

 allocate(Cin(kz%n,ksh%n,kxm%n))
 allocate(Cout(kz%n,ksh%n,kxm%n)); Cout=0;
 allocate(tfld2(ksh%n,kxm%n));tfld2=0;
 allocate(tfld1(kz%n));tfld1=0;
 allocate(temp(ksh%n));temp=0;
 allocate(krow(ksh%n));krow=0;
 allocate(shift(ksh%n));shift=0;

!--------------------------------------------------
! FFTW1D
 call fftw1d(tfld1)
 do ih=1,ksh%n
    do ix=1,kxm%n
       tfld1=cmplx(incig(:,ih,ix),0.)
       call fftexe1(tfld1)
       Cin(:,ih,ix)=tfld1/sqrt(1.*kz%n)
    end do
end do

!--------------------------------------------------
! FFTW2D
call fftw2d(tfld2)
 do iz=1,kz%n
    tfld2=Cin(iz,:,:)
    call fftexe2(tfld2)
    Cin(iz,:,:)=tfld2/sqrt(1.*ksh%n*kxm%n)
end do

   do ix=1,kxm%n
      do iz=1,kz%n

!         if (abs(kz%all(iz)).le.tiny(1.)) then
!            if (abs(kxm%all(ix)).le.tiny(1.)) then
!               ta=-kxm%d/kz%d*sign(1.,kz%all(iz))*sign(1.,kxm%all(ix))
!            else if (kxm%all(ix).lt.0.) then
!               ta=huge(1.)*sign(1.,kz%all(iz))
!            else if (kxm%all(ix).gt.0.) then
!               ta=-huge(1.)*sign(1.,kz%all(iz))
!            end if
!         else
!            ta=-kxm%all(ix)/kz%all(iz)
!         end if

         if (abs(kz%all(iz)).le.tiny(1.)) then
            ta=signal*sign(huge(1.),-kxm%all(ix))
         else
            ta=-kxm%all(ix)/kz%all(iz)
         end if


!----------------
        do ih=1,int(ksh%n/2)+1
          temp(ih+int(ksh%n/2))=Cin(iz,ih,ix)
        end do
        do ih=int(ksh%n/2)+2,ksh%n
          temp(ih-int(ksh%n/2)-1)=Cin(iz,ih,ix)
        end do

          shift=ksh%all*cig%sh%o
          temp=temp*cexp(cmplx(0,shift))
!----------------


         if (push) then

         do ih=1,ksh%n
  
            if (abs(kz%all(iz)).le.tiny(1.)) then
               khrot=signal*sign(ksh%all(ih),-kxm%all(ix))
            else
               tg=-ksh%all(ih)/kz%all(iz)
               if (abs(1-signal*ta*tg)>tiny(1.)) then
                  tteta=(ta+signal*tg)/(1.-signal*ta*tg)
                  khrot=-signal*kz%all(iz)*tteta
               else
                  khrot=signal*sign(ksh%all(ih),-kxm%all(ix))
               end if
            end if

               ikh=int(khrot/ksh%d+float(ksh%n/2))

               f=(khrot/ksh%d)+ksh%n/2 - ikh
               g=1-f
  
            if(ikh.ge.1.and.ikh.lt.ksh%n) then
               krow(ikh)= g*temp(ih) 
               krow(ikh)= f*temp(ih)
            else
               if(ikh.eq.ksh%n) then
                  krow(ikh)= g*temp(ih) 
                  krow(1)= f*temp(ih)
               else
                  if(ikh.eq.0) then
                     krow(ksh%n)=g*temp(ih) 
                     krow(1)=f*temp(ih)
                  end if
               end if
            end if
  
         end do

!----------------
          shift=-ksh%all*cig%sh%o
          krow=krow*cexp(cmplx(0,shift))

        do ih=1,int(ksh%n/2)+1
          Cout(iz,ih,ix)=krow(ih+int(ksh%n/2)-1)
        end do
        do ih=int(ksh%n/2)+2,ksh%n
          Cout(iz,ih,ix)=krow(ih-int(ksh%n/2)-1)
        end do
!----------------

!      end do
!   end do

   write(0,*) k1,k2

else

!   do ix=1,kxm%n
      do ih=1,ksh%n
!         do iz=1,kz%n
  
!            if (abs(kz%all(iz)).le.tiny(1.)) then
!               if (abs(kxm%all(ix)).le.tiny(1.)) then
!                  ta=-kxm%d/kz%d*sign(1.,kz%all(iz))*sign(1.,kxm%all(ix))
!               else if (kxm%all(ix).lt.0.) then
!                  ta=huge(1.)*sign(1.,kz%all(iz))
!               else if (kxm%all(ix).gt.0.) then
!                  ta=-huge(1.)*sign(1.,kz%all(iz))
!               else if (kxm%all(ix).eq.0.) then
!                  ta=0.
!               end if
!            else
!               ta=-kxm%all(ix)/kz%all(iz)
!            end if
  
            if (abs(kz%all(iz)).le.tiny(1.)) then
               if (abs(ksh%all(ih)).le.tiny(1.)) then
                  tg=-ksh%d/kz%d*sign(1.,kz%all(iz))*sign(1.,ksh%all(ih))
               else if (ksh%all(ih).lt.0.) then
                  tg=huge(1.)*sign(1.,kz%all(iz))
               else if (ksh%all(ih).gt.0.) then
                  tg=-huge(1.)*sign(1.,kz%all(iz))
               else if (ksh%all(ih).eq.0.) then
                  tg=0.
               end if
            else
               tg=-ksh%all(ih)/kz%all(iz)
            end if
  
            if (abs(1.-signal*ta*tg).le.tiny(1.)) then
               if (abs(tg+signal*ta).le.tiny(1.)) then
!                  tteta=huge(1.)
               else if (tg+signal*ta.lt.0) then
                  tteta=huge(1.)*sign(1.,1.-signal*ta*tg)
               else if (tg+signal*ta.gt.0) then
                  tteta=huge(1.)*sign(1.,1.-signal*ta*tg)
               else if (tg+signal*ta.eq.0) then
                  tteta=0.
               end if
            else
               tteta=(ta+signal*tg)/(1.-signal*ta*tg)
            end if
  
            khrot= -signal*kz%all(iz)*tteta
  
            if ((khrot.ge.ksh%all(1)).and.(khrot.le.ksh%all(int(ksh%n/2)+1))) then
               ikh=int((khrot-ksh%all(1))/ksh%d)+1.5
            else if ((khrot.ge.ksh%o).and.(khrot.le.ksh%all(ksh%n))) then
               ikh=int((khrot-ksh%o)/ksh%d)+1.5
            end if
  
            if (ikh.ge.1.and.ikh.le.ksh%n) then
               Cout(iz,ikh,ix)=Cout(iz,ikh,ix)+Cin(iz,ih,ix)
            end if
  
         end do
end if
      end do
   end do

!end if
!--------------------------------------------------
! IFFTW2D
 do iz=1,kz%n
    tfld2=Cout(iz,:,:)
    call ifftexe2(tfld2)
    Cout(iz,:,:)=tfld2/sqrt(1.*ksh%n*kxm%n)
end do

!--------------------------------------------------
! IFFTW1D

 do ih=1,ksh%n
    do ix=1,kxm%n
       tfld1=Cout(:,ih,ix)
       call ifftexe1(tfld1)
       outcig(:,ih,ix)=real(tfld1/sqrt(1.*kz%n))
    end do
 end do

end subroutine

end module


! allocate(tfld3(kz%n,ksh%n,kxm%n));tfld3=0;
!!--------------------------------------------------
!! FFTW3D
! call fftw3d(tfld3)
!write(0,*) "passei 01"
! tfld3=cmplx(incig,0.)
!write(0,*) "passei 02"
! call fftexe3(tfld3)
!write(0,*) "passei 03"
! Cin=tfld3/sqrt(1.*kz%n*ksh%n*kxm%n)
! tfld3=Cout
! call ifftexe3(tfld3)
! outcig=tfld3/sqrt(1.*kz%n*ksh%n*kxm%n)
