!====================================================================d
!   
!     GPE :  
!          Fast Oscillation in time t due to small epsilon
!     
!		i*u_t= H_xy + H_z/epsilon + beta |u|^2 u
!               u(t=0,x) = u_0(x), 0 < t < TFinal
!
!      Or equivally can be written in slow change variable "s = t/epsilon "
!
!		i*u_s = H_z + epsilon *(H_xy + beta |u|^2) u
!               u(s=0,x) = u_0(x), 0< s < TFinal/epsilon
!
!  Operators :  H_xy = -(\prt_xx+\prt_yy) ; H_z = -\prt_zz
!
!   File related description:
!       split.f90 contains 2nd and 4th order splitting method for slow changing  tau
!
!=====================================================================d
      program main
	use constant
	use fft_mkl
	implicit none 
      complex(8) :: psi(0:nx,0:ny)
      complex(8) :: psi_gd(0:nx,0:ny)
      complex(8) :: psit(0:nx,0:ny)
      complex(8) :: k1(0:nx,0:ny),k2(0:nx,0:ny)
      complex(8) :: k3(0:nx,0:ny),k4(0:nx,0:ny)

	real(8) :: en,en1,enkin,enpot,enint,enexc,entotal,enchem
	real(8) :: time

	integer(4) :: inds,i,j,k,kth,inloop,ind_total 

      call initial ! x,umx, vtemp 
	call initial_2d(psi)
	call initial_fft_mkl
     
      inds = 0
      time = 0.d0

	kth = 1
	ind_total = nint(TFinal/dtMacro)
      open(36,file='energy.dat',status='unknown')
	open(10,file='modex.bin',form='unformatted',access='direct',recl=32,Status='replace')
	open(20,file='modey.bin',form='unformatted',access='direct',recl=32,Status='replace')

	call  energy(enkin,enint,entotal,psi)
      write(36,40) time,entotal,enkin,enint	
      write(*,40) time,entotal,enkin,enint	

	psi_gd = psi
	call dfft_2d_mkl(psi_gd)
	write(10,rec=kth) psi_gd(0:7,0)
	write(20,rec=kth) psi_gd(0,0:7)
	kth = kth + 1

      do while(inds .lt.ind_total)
		inds = inds+1
		time = time + dtMacro
 
		!print *, inds,time
		call gradient(psi,psi_gd)
		k1 = dtMacro*psi_gd

		call  gradient(psi+0.5d0*k1,psi_gd)
		k2 = dtMacro*psi_gd

		call  gradient(psi+0.5d0*k2,psi_gd)
		k3 = dtMacro*psi_gd

		call  gradient(psi+k3,psi_gd)
		k4 = dtMacro*psi_gd

		psi = psi+ (k1+2.d0*k2+2.d0*k3+k4)/6.d0

		call  energy(enkin,enint,entotal,psi)
      	write(36,40) time,entotal,enkin,enint	
      	write(*,40) time,entotal,enkin,enint	
!
		psi_gd = psi
		call dfft_2d_mkl(psi_gd)
		write(10,rec=kth) psi_gd(0:7,0)
		write(20,rec=kth) psi_gd(0,0:7)
		kth = kth + 1
	enddo 
40	format(t4,4(2x,f16.10))

	close(10)
	close(20)

	call free_fft_mkl

101	format(t4,4(2x,f16.10))
	open(10,file='psi.dat',Status='unknown')
	do j =0,ny
		do i = 0,nx
			write(10,101) dreal(psi(i,j)),dimag(psi(i,j))
		enddo
	enddo
	close(10)

      write(*,*)  'it calculated for',inds,'steps'

       end 

!=====================================================================d
      subroutine  energy(enkin,enint,entotal,psi)
!======================================================================d
	use constant
	use fft_mkl
      implicit none
      complex(8),intent(in) :: psi(0:nx,0:ny)
	real(8),intent(out) :: enkin,enint,entotal

      complex(8) :: psit(0:nx,0:ny)

	integer(4) :: i,j
	real(8) :: t1,t3,hxy,hxyl

	psit = psi
	call dfft_2d_mkl(psit)
      hxyl = (xmax-xmin)*(ymax-ymin)
      hxy = hx*hy
      entotal = 0.0d0
      enkin = 0.d0
      enint = 0.0d0
      do 300 j=0,ny-1
            do 300 i=0,nx-1
       t1= abs(psit(i,j))**2.d0*vtemp(i,j)*hxyl
       t3= 0.5d0*beta*eps*abs(psi(i,j))**4.d0*hxy
       enkin = enkin + t1
       enint = enint + t3
       entotal = entotal + t1+t3
300	continue

      end
!=====================================================================d

