module fft_mkl
	use constant,only : nx,ny
	use mkl_dfti
	
	implicit none
	type(dfti_descriptor), POINTER :: dfft_x,dfft_y,dfft_z 
	integer(4) :: Status
	contains
		subroutine initial_fft_mkl
		Status = DftiCreateDescriptor(dfft_x, DFTI_DOUBLE,DFTI_COMPLEX,1,nx)
		status = DftiSetValue(dfft_x,DFTI_FORWARD_SCALE,1.d0/dble(nx)) 

		Status = DftiCreateDescriptor(dfft_y, DFTI_DOUBLE,DFTI_COMPLEX,1,ny)
		status = DftiSetValue(dfft_y,DFTI_FORWARD_SCALE,1.d0/dble(ny)) 

		Status = DftiCommitDescriptor(dfft_x)
		Status = DftiCommitDescriptor(dfft_y)

		end subroutine initial_fft_mkl
		subroutine free_fft_mkl

		Status = DftiFreeDescriptor(dfft_x)
		Status = DftiFreeDescriptor(dfft_y)

		end  subroutine free_fft_mkl

end module fft_mkl

	subroutine dfft_2d_mkl(psi)
!======================================================================d
!        inverse disrete sine transform in 3d
!======================================================================d
	use fft_mkl
	implicit none
	complex(8) :: psi(0:nx,0:ny)
	
	integer(4) :: i,j,k

	do j=0,ny-1
		Status = DftiComputeForward(dfft_x,psi(0:nx-1,j))
	enddo

	do i=0,nx-1
		Status = DftiComputeForward(dfft_y, psi(i,0:ny-1))
	enddo

	end subroutine dfft_2d_mkl

	subroutine difft_2d_mkl(psi)
	use fft_mkl
	implicit none
	complex(8) :: psi(0:nx,0:ny)
	integer(4) :: i,j

	do j=0,ny-1
		Status = DftiComputeBackward(dfft_x, psi(0:nx-1,j))
	enddo

	do i=0,nx-1
		Status = DftiComputeBackward(dfft_y, psi(i,0:ny-1))
	enddo

	end subroutine difft_2d_mkl


	subroutine gradient(psi,psi_gd)
	use constant, only : nx,ny,Per,eps,order
	implicit none
	complex(8),intent(in)  :: psi(0:nx,0:ny)
	complex(8),intent(out) :: psi_gd(0:nx,0:ny)

	complex(8) :: psi1(0:nx,0:ny)
	complex(8) :: psi2(0:nx,0:ny)
	complex(8) :: psi3(0:nx,0:ny)
	complex(8) :: psi4(0:nx,0:ny)
	integer(4) :: forward,backward

	integer(4) :: i,j

	forward = 1
	backward = 0
	if(order.eq.2) then 

		psi1 = psi
		psi2 = psi
		call split_2nd_2d(psi1,forward)
		call split_2nd_2d(psi2,backward)
		do j = 0,ny-1
			do i =0,nx-1
			psi_gd(i,j) = (psi1(i,j)-psi2(i,j))/(2.d0*Per*eps)
			enddo
		enddo

	elseif(order.eq.4) then

		do j = 0,ny-1
			do i =0,nx-1
			psi1(i,j) = psi(i,j)
			enddo
		enddo
		call split_4th_2d(psi1,forward)
		do j = 0,ny-1
			do i =0,nx-1
			psi2(i,j) = psi1(i,j)
			enddo
		enddo
		call split_4th_2d(psi2,forward)

		do j = 0,ny-1
			do i =0,nx-1
			psi_gd(i,j) = -psi2(i,j)+8.d0*psi1(i,j)
			enddo
		enddo

		do j = 0,ny-1
			do i =0,nx-1
			psi1(i,j) = psi(i,j)
			enddo
		enddo
		call split_4th_2d(psi1,backward)

		do j = 0,ny-1
			do i =0,nx-1
			psi2(i,j) = psi1(i,j)
			enddo
		enddo
		call split_4th_2d(psi2,backward)

		do j = 0,ny-1
			do i =0,nx-1
			psi_gd(i,j) = (psi_gd(i,j)-8.d0*psi1(i,j)+psi2(i,j))/(1.2d1*Per*eps)
			enddo
		enddo


	elseif(order.eq.6) then

		psi1 = psi
		call split_4th_2d(psi1,forward)
		psi2 = psi1
		call split_4th_2d(psi2,forward)
		psi3 = psi2
		call split_4th_2d(psi3,forward)
		psi_gd = psi3-9.d0*psi2+4.5d1*psi1

		psi1 = psi
		call split_4th_2d(psi1,backward)
		psi2 = psi1
		call split_4th_2d(psi2,backward)
		psi3 = psi2
		call split_4th_2d(psi3,backward)
		psi_gd = (psi_gd-psi3+9.d0*psi2-4.5d1*psi1)/(6.d1*Per*eps)

	elseif(order.eq.8) then

		psi1 = psi
		call split_2nd_2d(psi1,forward)
		psi2 = psi1
		call split_2nd_2d(psi2,forward)
		psi3 = psi2
		call split_2nd_2d(psi3,forward)
		psi4 = psi3
		call split_2nd_2d(psi4,forward)
		psi_gd = 3.d0*psi4-3.2d1*psi3+1.68d2*psi2-6.72d2*psi1

		psi1 = psi
		call split_2nd_2d(psi1,backward)
		psi2 = psi1
		call split_2nd_2d(psi2,backward)
		psi3 = psi2
		call split_2nd_2d(psi3,backward)
		psi4 = psi3
		call split_2nd_2d(psi4,backward)
		psi_gd = (psi_gd-(3.d0*psi4-3.2d1*psi3+1.68d2*psi2-6.72d2*psi1))/(-8.4d2*Per*eps)

	endif

	end subroutine

	subroutine split_2nd_2d(psi,flag)
	use constant,only: im, nx,ny,umx,umy,Per,dtMicro, Vtemp,beta,eps
	use fft_mkl
	implicit none
	complex(8) :: psi(0:nx,0:ny)
	integer(4) :: flag

	integer(4) :: i,j,k,ntMicro

	ntMicro = nint(Per/dtMicro)

	if(flag.eq.1) then
		do k = 1,ntMicro

			do j=0,ny-1
				do i = 0,nx-1
				psi(i,j) = exp(-im*beta*eps*dtMicro/2.d0*abs(psi(i,j))**2.d0)*psi(i,j)
				enddo
			enddo

		call dfft_2d_mkl(psi)
			do j=0,ny-1
				do i = 0,nx-1
				psi(i,j) =exp(-im*dtMicro*Vtemp(i,j))*psi(i,j)
				enddo
			enddo
		call difft_2d_mkl(psi)

			do j=0,ny-1
				do i = 0,nx-1
				psi(i,j) = exp(-im*beta*eps*dtMicro/2.d0*abs(psi(i,j))**2.d0)*psi(i,j)
				enddo
			enddo

		enddo
	else
		do k = 1,ntMicro
			do j=0,ny-1
				do i = 0,nx-1
				psi(i,j) = exp(im*beta*eps*dtMicro/2.d0*abs(psi(i,j))**2.d0)*psi(i,j)
				enddo
			enddo

		call dfft_2d_mkl(psi)
			do j=0,ny-1
				do i = 0,nx-1
				psi(i,j) = exp(im*dtMicro*Vtemp(i,j))*psi(i,j)
				enddo
			enddo
		call difft_2d_mkl(psi)

			do j=0,ny-1
				do i = 0,nx-1
				psi(i,j) = exp(im*beta*eps*dtMicro/2.d0*abs(psi(i,j))**2.d0)*psi(i,j)
				enddo
			enddo

		enddo

	endif

	end subroutine split_2nd_2d

	subroutine split_4th_2d(psi,flag)
	use constant
	use fft_mkl
	implicit none
	complex(8) :: psi(0:nx,0:ny)
	integer(4) :: flag

	integer(4) :: i,j,k,ntMicro
	ntMicro = nint(Per/dtMicro)

	if(flag.eq.1) then
		do k = 1,ntMicro
			do j=0,ny-1
				do i = 0,nx-1
			psi(i,j)= exp(-im*beta*eps*2.d0*w1*dtMicro*abs(psi(i,j))**2.d0)*psi(i,j)
				enddo
			enddo

		call dfft_2d_mkl(psi)
			do j=0,ny-1
				do i = 0,nx-1
			psi(i,j) = exp(-im*2.d0*w2*dtMicro*Vtemp(i,j))*psi(i,j)
				enddo
			enddo
		call difft_2d_mkl(psi)

			do j=0,ny-1
				do i = 0,nx-1
			psi(i,j) = exp(-im*beta*eps*2.d0*w3*dtMicro*abs(psi(i,j))**2.d0)*psi(i,j)
				enddo
			enddo

		call dfft_2d_mkl(psi)
			do j=0,ny-1
				do i = 0,nx-1
			psi(i,j) = exp(-im*2.d0*w4*dtMicro*Vtemp(i,j))*psi(i,j)
				enddo
			enddo
		call difft_2d_mkl(psi)

			do j=0,ny-1
				do i = 0,nx-1
			psi(i,j) = exp(-im*beta*eps*2.d0*w3*dtMicro*abs(psi(i,j))**2.d0)*psi(i,j)
				enddo
			enddo

		call dfft_2d_mkl(psi)
			do j=0,ny-1
				do i = 0,nx-1
			psi(i,j) = exp(-im*2.d0*w2*dtMicro*Vtemp(i,j))*psi(i,j)
				enddo
			enddo
		call difft_2d_mkl(psi)

			do j=0,ny-1
				do i = 0,nx-1
			psi(i,j) = exp(-im*beta*eps*2.d0*w1*dtMicro*abs(psi(i,j))**2.d0)*psi(i,j)
				enddo
			enddo

		enddo

	else
		do k = 1,ntMicro
			do j=0,ny-1
				do i = 0,nx-1
			psi(i,j)= exp(im*beta*eps*2.d0*w1*dtMicro*abs(psi(i,j))**2.d0)*psi(i,j)
				enddo
			enddo

		call dfft_2d_mkl(psi)
			do j=0,ny-1
				do i = 0,nx-1
			psi(i,j) = exp(im*2.d0*w2*dtMicro*Vtemp(i,j))*psi(i,j)
				enddo
			enddo
		call difft_2d_mkl(psi)

			do j=0,ny-1
				do i = 0,nx-1
			psi(i,j) = exp(im*beta*eps*2.d0*w3*dtMicro*abs(psi(i,j))**2.d0)*psi(i,j)
				enddo
			enddo

		call dfft_2d_mkl(psi)
			do j=0,ny-1
				do i = 0,nx-1
			psi(i,j) = exp(im*2.d0*w4*dtMicro*Vtemp(i,j))*psi(i,j)
				enddo
			enddo
		call difft_2d_mkl(psi)

			do j=0,ny-1
				do i = 0,nx-1
			psi(i,j) = exp(im*beta*eps*2.d0*w3*dtMicro*abs(psi(i,j))**2.d0)*psi(i,j)
				enddo
			enddo

		call dfft_2d_mkl(psi)
			do j=0,ny-1
				do i = 0,nx-1
			psi(i,j) = exp(im*2.d0*w2*dtMicro*Vtemp(i,j))*psi(i,j)
				enddo
			enddo
		call difft_2d_mkl(psi)

			do j=0,ny-1
				do i = 0,nx-1
			psi(i,j) = exp(im*beta*eps*2.d0*w1*dtMicro*abs(psi(i,j))**2.d0)*psi(i,j)
				enddo
			enddo
		enddo
	endif

	end subroutine split_4th_2d

