program gp
        use params
        use output
        !use utils
        !use rhs
        !use potential
        implicit none
        character(len=250) fname
        integer :: i, j, k, l, ii, jj, kk, ll
        double precision :: ret,n0,radius_r,rn1,rn2,norm1
        CALL init_params

        if (restart == 0) then 

          write(fname, '(a)') 'utils.0'
          open (8, FILE = fname)

          TIME = 0.0d0
          DT = -EYE*DTSIZE
                !$OMP PARALLEL DO
                do i = -NX/2,NX/2
                do j = -NY/2,NY/2
                do k = -NZ/2,NZ/2
                do l = -NW/2,NW/2
                                radius_r =sqrt((dble(i)*DSPACE)**2.0d0+(dble(j)*DSPACE)**2.0d0&
                                        +(dble(k)*DSPACE)**2.0d0+(dble(l)*DSPACE)**2.0d0)
                                if (radius_r.lt.POT_RADIUS) then
                                        OBJPOT(i,j,k,l) = 0.0d0
                                else
                                        OBJPOT(i,j,k,l) = 10.0d0
                                end if
                end do
                end do
                end do
                end do
                !$OMP END PARALLEL DO

                !$OMP PARALLEL DO
                do i = -NX/2, NX/2
                do j = -NY/2, NY/2
                do k = -NZ/2, NZ/2
                do l = -NW/2, NW/2
                                        GRID(i,j,k,l)=(1.0d0,0.0d0)-REAL(OBJPOT(i,j,k,l))
                                        if (REAL(OBJPOT(i,j,k,l)).gt.1.0d0) then
                                                GRID(i,j,k,l) = 0.0d0 +EYE*0.0d0
                                        end if
                end do
                end do
                end do
                end do
                !$OMP END PARALLEL DO
                call int_grid(CONJG(GRID)*GRID,NORM)
         
          call runit_im(0,0)
          include 'ic.in'
          call runit_v_toll(0,0)
          
          DT = DTSIZE
          
          TIME = 0.0d0
          
          if (rot_real) then
            omega_rot=omega_rot_initial
          end if
          
          
          TIME = 0.0d0
          !call add_noise
          
          call runit(NSTEPS,1,1)
          
          close(8)
          
          close(9)
        
        end if

        if (restart == 1) then

                DT = DTSIZE
                TIME = DTSIZE*dble(OLD_STEPS)

                write(fname, '(a,i0.4)') '0000.dumpwf.',IMAGSTART/DUMPWF
                call read_wf_file_nopot(fname)
                write(6,"(a)") fname 
                !call calc_misc(energy_gs)
                !$OMP PARALLEL DO
                do i = -NX/2,NX/2
                do j = -NY/2,NY/2
                do k = -NZ/2,NZ/2
                do l = -NW/2,NW/2
                                radius_r=sqrt((dble(i)*DSPACE)**2.0d0+(dble(j)*DSPACE)**2.0d0&
                                        +(dble(k)*DSPACE)**2.0d0+(dble(l)*DSPACE)**2.0d0)
                                if (radius_r.lt.POT_RADIUS) then
                                        OBJPOT(i,j,k,l) = 0.0d0
                                else
                                        OBJPOT(i,j,k,l) = 10.0d0
                                end if
                end do
                end do
                end do
                end do
                !$OMP END PARALLEL DO   
                call runit(NSTEPS,1,1)
        
                close(8)
                close(9)              
         end if

        if (restart == 2) then 

                write(fname, '(a)') 'utils.0'
                open(8, FILE = fname)
                open(9, FILE='norm.dat')
          
                TIME = 0.0d0
  

                DT = -EYE*DTSIZE
                call calc_OBJPOT 
                write(fname, '(a)') '0000.dumpwf_ground_state.0000'
                call read_wf_file_nopot(fname)  
                call dump_wavefunction(0,0)        
                !call calc_misc(energy_gs)
                !$OMP PARALLEL DO
                do i = -NX/2,NX/2
                do j = -NY/2,NY/2
                do k = -NZ/2,NZ/2
                do l = -NW/2,NW/2
                                radius_r=sqrt((dble(i)*DSPACE)**2.0d0+(dble(j)*DSPACE)**2.0d0&
                                        +(dble(k)*DSPACE)**2.0d0+(dble(l)*DSPACE)**2.0d0)
                                if (radius_r.lt.POT_RADIUS) then
                                        OBJPOT(i,j,k,l) = 0.0d0
                                else
                                        OBJPOT(i,j,k,l) = 10.0d0
                                end if
                end do
                end do
                end do
                end do
                !$OMP END PARALLEL DO  
          call int_grid(CONJG(GRID)*GRID,NORM)
          write(6,*) "norm before = ", NORM
          include 'ic.in'
          call int_grid(CONJG(GRID)*GRID,NORM1)
          write(6,*) "norm after = ",NORM1
          GRID = sqrt(NORM)*GRID/sqrt(NORM1)
          call dump_wavefunction(0,5)
          call runit_v_toll(0,0)
          DT = DTSIZE
          TIME = 0.0d0

          if (rot_real) then
	    omega_rot=omega_rot_initial
	  end if
	  
	  TIME = 0.0d0
	  call runit(NSTEPS,1,1)
	  
	  close(8)
          close(9)

        end if

end PROGRAM gp

subroutine runit(steps,rt,plot)
        use params
        implicit none
        integer :: steps,rt,i,plot, remeshcount,ii,jj,kk,ll,ramsey_pause
        double precision :: energy,energyi,NORM1
        double precision :: rx,ry,rz,rn1,rn2
        write (unit=8,fmt="(a)") 'real time advancement'
        flush(8)

                        call dump_wavefunction(0,0) !dump initial condition

        do i = OLD_STEPS+1, steps
                call iterate(rt)
                TIME = TIME + dble(DT)
                if (modulo(i,dumputil) == 0) then
                        if (plot == 1) then
                                !call calc_misc(energy)
                                !call calc_norm
                                call int_grid(CONJG(GRID)*GRID,NORM)
                                write (unit=6,fmt="(a,f15.8,a,i5,a,f15.6)")&
                                "Time : ",TIME, " iteration ",i,", N:  ",NORM
                                flush(6)
                        end if
                end if

                if (modulo(i,10) == 0) then
                        open (10, FILE = "STATUS")
                        if (rt == 1) then
                                write (unit=10,fmt="(a,f6.2,a)")&
                                "Simulating::     ",(dble(i)/dble(steps))*100.0d0,"%"
                                else
                                        write (unit=10,fmt="(a,f6.2,a)")&
                                        "Ground State:     " ,(dble(i)/dble(steps))*100.0d0,"%"
                        end if
                close(10)
                end if
                if (modulo(i,dumpwf) == 0) then
                        if (rt == 1) then
                                if (plot == 1) then
                                        call dump_wavefunction(i,0)
                                end if
                        end if
                end if

        end do
end subroutine

subroutine runit_im(rt,plot)
	use params
	implicit none
	integer :: rt,i,plot,j,k,l
        double precision :: energy_old,energy,err
        i=0
        write(6,*) i
        !energy_old=10.0d0
        energy=1.0d0
!        call calc_OBJPOT
!        call initCond

	write (unit=8,fmt="(a)") 'imaginary time advancement'
	flush(8)

	write (unit=9,fmt="(a)") 'imaginary time advancement'
	flush(9)
	call calc_misc(energy)
	energy_old=10000.0d0!*!energy
        do i = 1,100!while ((abs(energy-energy_old)/energy_old).gt.(toll))
                !i=i+1
                energy = 1.0d0!(abs(energy-energy_old)/energy_old)
                call iterate(rt)
                TIME = TIME + DIMAG(DT)
                if(MODULO(i,dumputil)==0) then
                write(6,*)   "Time : ",TIME, " iteration ",i,", mid density:",real(conjg(GRID(0,0,0,0))*(GRID(0,0,0,0)))
                end if
                if (modulo(i,dumpwf_imag) == 0) then
                        call dump_wavefunction(i,3)
                end if
                energy_old=energy
                flush(6)
                call calc_misc(energy)
        end do
        energy_gs=energy
	call dump_wavefunction(0,1)
	
end subroutine

subroutine runit_v(steps,rt,plot)
	use params
	implicit none
	integer :: steps,rt,i,plot
	double precision :: energy

	write (unit=8,fmt="(a)") 'imaginary time advancement with vortices'
	flush(8)

	write (unit=9,fmt="(a)") 'imaginary time advancement with vortices'
	flush(9)

	do i = 1, steps
		call iterate(rt)
		TIME = TIME + dble(DT)
		if (modulo(i,dumputil_vortex) == 0) then
		!	call calc_misc(energy)
		end if
		if (modulo(i,10) == 0) then
			open (10, FILE = "STATUS")
			if (rt == 1) then
					write (unit=10,fmt="(a,i3,a,f6.2,a)")&
						"Simulating: ",VOBS, ":     ",(dble(i)/dble(steps))*100.0d0,"%"
				else
					write (unit=10,fmt="(a,i3,a,f6.2,a)")&
						"Ground State: ",VOBS, ":     " ,(dble(i)/dble(steps))*100.0d0,"%"
			end if
			close(10)
		end if
		if (modulo(i,dumpwf_vortex) == 0) then
			call dump_wavefunction(i,2)
		end if
		if(potRep .eq. 1 .and. rt == 1) then
			call calc_OBJPOT
		end if
	end do
end subroutine


subroutine runit_back(rt,plot)
        use params
        implicit none
        integer :: rt,i,plot
        double precision :: energy, energy_old,err

        i=0
        call dump_wavefunction(i,2)
        energy_old=energy_gs

        write (unit=6,fmt="(a)") 'imaginary time advancement with vortices'
        flush(6)
        DT = -1.0d0*DTSIZE
        TIME = 0.0d0
        call calc_misc(energy)
        do i = 1,NSTEPS!while ((abs(energy-energy_old)/energy_old).gt.(toll_v))

!               i=i+1
!
                call iterate(rt)
                TIME = TIME + DT
                if (modulo(i,5) == 0) then
                                write (unit=6,fmt="(a,f7.3)") "Time:",TIME
                end if

                if (modulo(i,dumpwf_vortex) == 0) then
                        call dump_wavefunction(i,3)
                end if

        end do

end subroutine

subroutine runit_v_toll(rt,plot)
	use params
	implicit none
	integer :: rt,i,plot
	double precision :: energy, energy_old,err

	i=0
        call dump_wavefunction(i,2)
        energy_old=energy_gs

	write (unit=6,fmt="(a)") 'imaginary time advancement with vortices'
	flush(6)


        call calc_misc(energy)

	do i = 0,10!while ((abs(energy-energy_old)/energy_old).gt.(toll_v))

!		i=i+1
!
		call iterate(rt)
		TIME = TIME + DIMAG(DT)
                !norm
		     !     call calc_norm
        !NORM_1=NORM
        !GRID = sqrt(NORM_0)*GRID/sqrt(NORM_1)
        !norm
	!	energy_old=energy
		
	!	call calc_misc(energy)
         !       err = (abs(energy-energy_old)/energy_old)
         !       call calc_norm
         !       if(modulo(i,dumpwf_vortex)==0) then
         !       write (unit=13,fmt="(f8.4,3x,I0.3,3x,f15.12,3x,f15.6)")time,i/dumpwf_vortex, err, norm
         !       flush(13)
         !       end if
		if (modulo(i,5) == 0) then


				write (unit=6,fmt="(a,f7.3)") "Immaginary time:" ,TIME
			

		end if

		if (modulo(i,dumpwf_vortex) == 0) then
			call dump_wavefunction(i,3)
		end if

	end do



end subroutine

subroutine dump_wavefunction (II,typ)
	use params
	use output
	implicit none
	integer :: II,i,j,k,typ
	character(len=250) fname

	if (typ == 0) then
		write(fname, '(i0.4,a,i0.4)') VOBS,'.dumpwf.',II/dumpwf
	end if

	if (typ == 1) then
		write(fname, '(i0.4,a,i0.4)') VOBS,'.dumpwf_ground_state.',II/dumpwf
	end if

	if (typ == 2) then
		write(fname, '(i0.4,a,i0.4)') VOBS,'.dumpwf_vortex.',II/dumpwf_vortex
	end if
	
	if (typ == 3) then
		write(fname, '(i0.4,a,i0.4)') VOBS,'.dumpwf_imag.',II/dumpwf_imag
	end if
	
	if (typ == 4) then
		write(fname, '(i0.4,a,i0.4)') VOBS,'.dumpwf_diss.',II/dumpwf_diss
	end if
        if (typ == 5) then
                write(fname, '(i0.4,a,i0.4)') VOBS,'.dumpwf_add_vortices.',II/dumpwf
        end if
	call make_file_nopot(fname)
	call write_wf_file_nopot
end subroutine