PROGRAM GLOBAL
c
c	This program treats the evolution of an NDIM dimensional
c 	according to the non-linear sigma model. 
c	It calls subroutines that  evolves a background density field
c	(enersigma) and subroutines that evolve the sigma field 
c
c
c	This version of the program is VECTORIZED for the CONVEX
c
c  PARAMETERS:
c
c	The following parameters can be varied to study different
c	topological defects under different initial conditions:
c
c		NP,MP= determines the grid size.   
c		NDIM =	determines the number of internal dimensions
c			associated with the psi field
c		TX = determines the epoch of the transition from
c		      matter to radiation domination.
c		      the program assumes that the scale factor
c		      scales as a(t) = (t/tx) + 0.25*(t/tx)**2
c		      teq = the time of equal matter and radiation 
c			  =  0.828 tx
c		T0 = the starting time of the simulation
c		TSTOP = stoping time of the simulation
c		PRESSURE = determines whether code should compute the
c			    equation of state for the Goldstone bosons
c						(should be .false.)
c		NSTEP = maximum number of steps
c		IDUM = this is the seed for generating initial conditions
c			varying this parameter will generate a new realization
c			of the model
c
c  UNITS:
c	T:	time is measured in conformal time
c	R:	TX DETERMINES THE SIZE OF THE BOX:
c			the box size is np*(19.38 h^{-2} Mpc)/tx
c	PERT:   density fluctuations in the dark matter
c		 are normalized	in units of 1.5 \epsilon a(t_0)
c		 where \epsilon = 8\pi^2 G\phi_0^2  and a(t_0) is the
c		 scale factor at time t_0 and t_0 >> t_eq.
c
c
c INTERNAL VARIABLES:
c	
c	PSI(N,I,J,K) = the value of the Nth component of the field at i,j,k
c	PSIDOT(N,I,J,K) = the value of the time derivative
c	PERT(I,J,K) = the amplitude of the CDM perturbation at i,j,k
c	TAU00(I, J,K) = the amplitude of the pseudoenergy density at (i,j,k)
c	ENTROPY(I,J,K) = the amplitude of the entropy fluctuations
c	ENTROPYDOT(I,J,K) = its time derivative
c	NOLD =	number of old knots
c	NNEW = number of new knots
c	
C
C
C OUTPUT FILES:
C		xxx.scaling: records the scaling of the goldstone bosons
C		xxx.knots: records the number of knots
C		xxx.list: list the location of the knots 
C		xxx.pert.dat: density field (unformatted)
c
c
c		
c     Choose the simulation box size here : ndim =4 is
c     for texture, ndim = 3 is for global monopoles etc.
c     numerical grid is np times np times mp
c
	parameter (np=256,mp=256,ndim=4,pi=3.1415927,nfreq=100)
c
c	COMMON BLOCKS:
c
	dimension psi(np,ndim,np,mp),psid(np,ndim,np,mp)
	dimension fieldk(np,np,mp),fieldt(np,np,mp),fieldg(np,np,mp)
	dimension favk(np/2,np/2,mp/2)
	dimension favt(np/2,np/2,mp/2)
	dimension favg(np/2,np/2,mp/2)
	common w0,t0,tx,alpha,nu,beta
C
C	ARRAYS:
C
	dimension knotold(3,np*np),knotnew(3,np*np)
	dimension psihalf(np,ndim,np,mp),psipr(np,ndim,np,mp)
	dimension dpsidx(ndim),dpsidy(ndim),dpsidz(ndim),psiav(ndim)
	logical init
c    w0 a parameter which determines time from which density
c    perturbations are integrated. Normally take w0=0.
	w0=0.
	dt = 0.2
c starting simulations with correct scaling energy density
	if(ndim.eq.2) then
		t0=7.
	else if(ndim.eq.3) then
		t0=3.2
	else if(ndim.eq.4) then
		t0=2.
	else if(ndim.eq.6) then
		t0=1.
	endif
	t=t0
C
C Adjust tx to choose the size of simulation you need:
C tx= 19.38  corresponds to np h^-2 Mpc for example.
C
C box = size of simulation volume in h^-2 Mpc
C
c Caution: must run well into matter era to correctly
c compute coefficient of growing mode perturbation
c constant needed to normalise perturbations
c with epsilon = 10^-4 from COBE
c
	nstep = 100000
C seed for random number generator
	idum=-6457
C
C	OUTPUT FILES:
C
C
C	RECORD RUN INFORMATION
C
C	SET-UP THE INITIAL CONDITIONS:
C		VELINIT- gives the array a random initial value
c			 on the N-sphere.  (radius of N-sphere is 1)
C
C
	iunit=20
	call velinit(psi,ndim,np,mp,1.0,idum)
	do k = 1,mp
		do j=1,np
			do n = 1,ndim
				do i=1,np
				psipr(i,n,j,k)=psi(i,n,j,k)
				end do
			end do
		end do
	end do
	nold=0
	nwrite=2.01/dt
	write(6,*) nwrite
	flipflag=1
	do jstep=1,nstep
C	advance the psi field by one grid unit
C
		if(t.lt.3.*np/2.) call sstep(np,mp,ndim,psi,psipr
     >			,psid,psihalf,t,dt,flipflag)
c invoking `spin flip' to unwind textures correctly for ndim =4
c		if(t.gt.5..and.ndim.eq.4) flipflag=-1
		write(*,*) ' Taken one step, t=',t
c

c
c	increment time 
c
		t=t+dt
c
c	exit program
c
c		if(mod(jstep,nwrite).eq.0) then
c	if((t.gt.50.).and.(mod(jstep,nwrite).eq.0)) then
	if(t.gt.50.) then
	do k = 1,mp
		do j=1,np
			do i=1,np
	fieldk(i,j,k)=0.0
	fieldt(i,j,k)=0.0
	fieldg(i,j,k)=0.0
			enddo
		end do
	end do
	do k = 1,mp
		do j=1,np
			do n=1,ndim
			do i=1,np
	fieldk(i,j,k)=fieldk(i,j,k)+(psipr(i,n,j,k)-psi(i,n,j,k))**2/dt/dt
			enddo
			enddo
		end do
	end do
	toptot=0.
	do k = 1,mp
		do j=1,np
			do i=1,np
         kp=k+1
         if(kp.gt.mp)kp=1
         km=k-1
         if(km.eq.0)km=mp
         jp=j+1
         if(jp.gt.mp)jp=1
         jm=j-1
         if(jm.eq.0)jm=mp
         ip=i+1
         if(ip.gt.mp)ip=1
         im=i-1
         if(im.eq.0)im=mp
  		topdens=0.
			do n=1,ndim
	fieldg(i,j,k)=fieldg(i,j,k) +(
     >     	(psi(ip,n,j,k)-psi(im,n,j,k))**2+
     >	(psi(i,n,jp,k)-psi(i,n,jm,k))**2+
     >	(psi(i,n,j,kp)-psi(i,n,j,km))**2)/4.
	  dpsidx(n)=( 
     >	 psi(ip,n,j,k)+psi(ip,n,jp,k)+psi(ip,n,j,kp)+psi(ip,n,jp,kp)
     >   -(psi(i,n,j,k)+psi(i,n,jp,k)+psi(i,n,j,kp)+psi(i,n,jp,kp))
     >          )/4.0
	  dpsidy(n)=( 
     >	 psi(i,n,jp,k)+psi(ip,n,jp,k)+psi(ip,n,jp,kp)+psi(i,n,jp,kp)
     >   -(psi(i,n,j,k)+psi(ip,n,j,k)+psi(i,n,j,kp)+psi(ip,n,j,kp))
     >          )/4.0
	  dpsidz(n)=( 
     >	 psi(i,n,j,kp)+psi(ip,n,j,kp)+psi(i,n,jp,kp)+psi(ip,n,jp,kp)
     >   -(psi(i,n,j,k)+psi(i,n,jp,k)+psi(ip,n,j,k)+psi(ip,n,jp,k))
     >          )/4.0
	     psiav(n)=psi(i,n,jp,k)+psi(ip,n,jp,k)+psi(ip,n,jp,kp)
     >       +psi(i,n,jp,kp)+
     >       psi(i,n,j,k)+psi(ip,n,j,k)+psi(i,n,j,kp)
     >       +psi(ip,n,j,kp)
	psiav(n)=psiav(n)/8.	
			enddo
c   calc top density at i,j,k (centre of cube)
	in1=1
	in2=2
	in3=3
	in4=4
	topdens=topdens
     >      +psiav(in1)*dpsidx(in2)*dpsidy(in3)*dpsidz(in4)
     >      -psiav(in1)*dpsidx(in2)*dpsidy(in4)*dpsidz(in3)
     >      -psiav(in1)*dpsidx(in3)*dpsidy(in2)*dpsidz(in4)
     >      +psiav(in1)*dpsidx(in3)*dpsidy(in4)*dpsidz(in2)
     >      +psiav(in1)*dpsidx(in4)*dpsidy(in2)*dpsidz(in3)
     >      -psiav(in1)*dpsidx(in4)*dpsidy(in3)*dpsidz(in2)
	in1=4
	in2=1
	in3=2
	in4=3
	topdens=topdens
     >     -(psiav(in1)*dpsidx(in2)*dpsidy(in3)*dpsidz(in4)
     >      -psiav(in1)*dpsidx(in2)*dpsidy(in4)*dpsidz(in3)
     >      -psiav(in1)*dpsidx(in3)*dpsidy(in2)*dpsidz(in4)
     >      +psiav(in1)*dpsidx(in3)*dpsidy(in4)*dpsidz(in2)
     >      +psiav(in1)*dpsidx(in4)*dpsidy(in2)*dpsidz(in3)
     >      -psiav(in1)*dpsidx(in4)*dpsidy(in3)*dpsidz(in2))
	in1=2
	in2=3
	in3=4
	in4=1
	topdens=topdens
     >     -(psiav(in1)*dpsidx(in2)*dpsidy(in3)*dpsidz(in4)
     >      -psiav(in1)*dpsidx(in2)*dpsidy(in4)*dpsidz(in3)
     >      -psiav(in1)*dpsidx(in3)*dpsidy(in2)*dpsidz(in4)
     >      +psiav(in1)*dpsidx(in3)*dpsidy(in4)*dpsidz(in2)
     >      +psiav(in1)*dpsidx(in4)*dpsidy(in2)*dpsidz(in3)
     >      -psiav(in1)*dpsidx(in4)*dpsidy(in3)*dpsidz(in2))
	in1=3
	in2=4
	in3=1
	in4=2
	topdens=topdens
     >      +psiav(in1)*dpsidx(in2)*dpsidy(in3)*dpsidz(in4)
     >      -psiav(in1)*dpsidx(in2)*dpsidy(in4)*dpsidz(in3)
     >      -psiav(in1)*dpsidx(in3)*dpsidy(in2)*dpsidz(in4)
     >      +psiav(in1)*dpsidx(in3)*dpsidy(in4)*dpsidz(in2)
     >      +psiav(in1)*dpsidx(in4)*dpsidy(in2)*dpsidz(in3)
     >      -psiav(in1)*dpsidx(in4)*dpsidy(in3)*dpsidz(in2)
	fieldt(i,j,k)=topdens/2./pi**2
	toptot=toptot+topdens/2./pi**2
			enddo
		end do
	end do
	write(15,*) t,toptot
		if(mod(jstep,nwrite).eq.-500) then
	do l=1,mp/2
		do m=1,np/2
			do n=1,np/2
	favk(l,m,n)=(fieldk(2*l-1,2*m-1,2*n-1)
     >  	    +fieldk(2*l-1+1,2*m-1,2*n-1)
     >  	    +fieldk(2*l-1,2*m-1+1,2*n-1)
     >  	    +fieldk(2*l-1,2*m-1,2*n-1+1)
     >  	    +fieldk(2*l-1+1,2*m-1+1,2*n-1)
     >  	    +fieldk(2*l-1+1,2*m-1,2*n-1+1)
     >  	    +fieldk(2*l-1,2*m-1+1,2*n-1+1)
     >  	    +fieldk(2*l-1+1,2*m-1+1,2*n-1+1))/8.
	favt(l,m,n)=(fieldt(2*l-1,2*m-1,2*n-1)
     >  	    +fieldt(2*l-1+1,2*m-1,2*n-1)
     >  	    +fieldt(2*l-1,2*m-1+1,2*n-1)
     >  	    +fieldt(2*l-1,2*m-1,2*n-1+1)
     >  	    +fieldt(2*l-1+1,2*m-1+1,2*n-1)
     >  	    +fieldt(2*l-1+1,2*m-1,2*n-1+1)
     >  	    +fieldt(2*l-1,2*m-1+1,2*n-1+1)
     >  	    +fieldt(2*l-1+1,2*m-1+1,2*n-1+1))/8.
	favg(l,m,n)=(fieldg(2*l-1,2*m-1,2*n-1)
     >  	    +fieldg(2*l-1+1,2*m-1,2*n-1)
     >  	    +fieldg(2*l-1,2*m-1+1,2*n-1)
     >  	    +fieldg(2*l-1,2*m-1,2*n-1+1)
     >  	    +fieldg(2*l-1+1,2*m-1+1,2*n-1)
     >  	    +fieldg(2*l-1+1,2*m-1,2*n-1+1)
     >  	    +fieldg(2*l-1,2*m-1+1,2*n-1+1)
     >  	    +fieldg(2*l-1+1,2*m-1+1,2*n-1+1))/8.
			enddo
		enddo
	enddo
	iunit=iunit+1
			open(unit,form='unformatted',
     > 				status='unknown')
			write(iunit) 
     &                  (((favk(i,j,k),i=1,np/2),j=1,np/2),k=1,mp/2)
			write(iunit) 
     &                  (((favt(i,j,k),i=1,np/2),j=1,np/2),k=1,mp/2)
			write(iunit) 
     &                  (((favg(i,j,k),i=1,np/2),j=1,np/2),k=1,mp/2)
			close(iunit)
			endif
		endif
	if(t.gt.150.) stop
	enddo
	stop
	end
C*****************************************************************************
C
C 			SUBROUTINE SSTEP
C
C written by Neil Turok, Spring 1992
C modified June 24, 1992 by Ue-Li Pen
C
C*****************************************************************************
c     evolve a scalar field according to the nonlinear sigma model.
C
c ARGUMENTS:
c(I)   integer np,np,mp:  x,y,z dimensions of psi
c(I)   real    psi   :  ndim*ngs size array holding the field psi
c(I/O) real    psipr :  psi at the previous timestep, which also becomes psi
c			at the next time step
c(I)   real    psid2(1,1,1): if this first element is 9876543, don't bother
c                          updating psid2.  (If it were C, I would test against
c                          NULL).
c(I)   real    ndim  :  # of components of the scalar field.  =4 for textures.
c(I)   real    dt    :  best at <= .2
c(I)   real    a1,a2 :  the expansion factors at the previous and next half
c                       steps, with some dt**2's.
c(I)   real    flipflag:  = 1 if you DON'T want to flip spins, and
c 	  		   -1 if you DO want to flip
c
c 1. key:  (I) = input only (not modified)
c          (O) = output only (not read)
c
c 2. note that t and dt are not changed
c     
c
c 4. usage note:  flipflag should be +1 for the first few steps until field
c              has settled.
c
c     
c 5. vectorizes completely.
c
c 6. Attn:  it is important that the subroutine be completely self-contained,
c        i.e., does not call anything nor use any common blocks outside
c        of this file.  It should communicate to the outside world only as
c       specified in the comments above.

      subroutine sstep(np,mp,ndim,psi,psipr,psid,psihalf,t,dt,flipflag)
      dimension psi(np,ndim,np,mp),psid(np,ndim,np,mp)
     &      ,psihalf(np,ndim,np,mp)
     &     ,psipr(np,ndim,np,mp)
      logical parity
      save parity
      data parity/.true./
       a1 = (aexp(t-dt/2.)/aexp(t+dt/2.))**2
       a2 =(aexp(t)/aexp(t+dt/2.))**2*dt*dt
  
      if (parity .eqv. .true.) then
         parity = .false.
	call sstep1(np,mp,ndim,psi,psipr,psid,psihalf,dt,a1,a2,flipflag)
      else
         parity = .true.
	call 
     &    sstep1(np,mp,ndim,psipr,psi,psid,psihalf,dt,a1,a2,flipflag)
      endif
      return
      end
C*******************************************************************************
c
c	This subroutine evolves the psi field by a modified leapfrog
c	algorithm
c
c
c June 24, 1992:
cupen:  I changed the arguments and vectorization arrays.
 
	subroutine sstep1(np,mp,ndim,
     &     psi,psipr,psid,psihalf,dt,a1,a2,flipflag)
      parameter (pi=3.14159265, MAXNDIM=12)
      dimension psi(np,ndim,np,mp),psid(np,ndim,np,mp)
     >   ,psihalf(np,ndim,np,mp)
     >    ,psipr(np,ndim,np,mp)
      dimension dpsi2(1024),psdpsi(1024),grc(1024),xroot(1024)
      dimension psitem(1024),ddpsi(1024),ekinp(1024),egrdp(1024)
      dimension dpsi(1024,MAXNDIM)
c     
c
c
c
c	these terms are used to advance psipr
c
c	a1 = ((t-dt/2.)**2/(t+dt/2.)**2)**2
c	a2 =((t)**2/(t+dt/2.)**2)**2*dt*dt
c	loop over grid
c
      dtinv = 1./dt
      if (np .gt. 1024) then
	write(*,*) "error: np>1024"
	stop
	endif
      if (ndim .gt. MAXNDIM) write(*,*) "error: ndim>MAXNDIM"
      one=1.
      do k=1,mp
         kp=k+1
         if(kp.gt.mp)kp=1
         km=k-1
         if(km.eq.0)km=mp
         do j=1,np
	    jp=j+1
	    if(jp.gt.np)jp=1
	    jm=j-1
	    if(jm.eq.0)jm=np
            do i=1,np
               dpsi2(i)=0.
               psdpsi(i)=0.
               grc(i)=0.
            enddo
	    do i=1,np
        	ip=i+1
            if(ip.gt.np)ip=1
            im=i-1
            if(im.eq.0)im=np
               do n=1,ndim
                  ddpsi(i)=psi(ip,n,j,k)+psi(im,n,j,k)+psi(i,n,jp,k)
     >                 +psi(i,n,jm,k)+psi(i,n,j,kp)+psi(i,n,j,km)
                  grc(i)=grc(i)+psi(i,n,j,k)*ddpsi(i)
                  dpsi(i,n)=psi(i,n,j,k)-a1*psipr(i,n,j,k)+a2*ddpsi(i)
                  dpsi2(i)=dpsi2(i)+ dpsi(i,n)**2
                  psdpsi(i)= psdpsi(i)+ psi(i,n,j,k)*dpsi(i,n)
               enddo
c               grc(i)=max(flipflag,sign(one,grc(i)))
          xroot(i)= sqrt(1.-dpsi2(i)+psdpsi(i)**2) - psdpsi(i)
            enddo
            do nu=1,ndim
                do i=1,np
               psipr(i,nu,j,k)= xroot(i)*psi(i,nu,j,k)+dpsi(i,nu)
	     psid(i,nu,j,k) = (psipr(i,nu,j,k)-psi(i,nu,j,k))
     >					*dtinv
	     psihalf(i,nu,j,k)=0.5*(psipr(i,nu,j,k)+psi(i,nu,j,k))
                  enddo
c
            enddo
         enddo
      enddo
c
c
      return
      end


	subroutine velinit(arr,ndim,np,mp,v,idum)
	real arr(np,ndim,np,mp),x(10)
	do k=1,mp
		do i=1,np
			do j=1,np
333				continue
				sum = 1.e-35
				do n=1,ndim
					x(n)=2.*ran1(idum)-1.
					sum = sum + x(n)*x(n)
				end do
				if(sum.gt.1.0) goto 333
				sum = sqrt(sum)
				vsum = v/sum
				do n=1,ndim
					arr(i,n,j,k) = x(n)*vsum
				end do
			enddo
		enddo
	enddo
	return
	end
C**************************************************************************
C
C	FUNCTION RAN1
C		Numerical Recipes Random Number Generator
C
C
C****************************************************************************
	function ran1(idum)
	integer idum,ia,im,iq,ir,ntab
	real ran1,am,atab	
	parameter (ia=16807,im=2147483647,am=1./im,iq=127773,ir=2836,
     *		ntab=32,atab=ntab-1)
c	integer j,k
	real y,v(ntab)
	data v /ntab*0./, y /0.5/
	if (idum.le.0) then
c
		idum=max(-idum,1)
c
		do j=ntab,1,-1
c
			k=idum/iq
			idum=ia*(idum-k*iq)-ir*k
			if (idum.lt.0) idum=idum+im
			v(j)=am*idum
		enddo
		y=v(1)
	endif
1	continue
c
		k=idum/iq
		idum=ia*(idum-k*iq)-ir*k
c
		if (idum.lt.0) idum=idum+im
		j=1+int(atab*y)
c
		y=v(j)
		ran1=y
c
		v(j)=am*idum
	if (ran1.eq.0..or.ran1.eq.1.)goto 1
c
	return
	end
c
c       FUNCTION AEXP(T)
c               computes the expansion factor
c
c
        function aexp(t)
        common w0,t0,tx,alpha,nu,beta
c        aexp = t/tx *(1.+ 0.25*(t/tx))
        aexp = t**4/(1+(t/5.)**2)
        return
        end