!8/22/01
! lame6.hpf
! STATUS: incomplete, workable
!  -must find cause of "segmentation fault"
!  -re-modifying the correlation to be CORRECTLY periodic- corrbf2 
!     and corrln have to be completed

! QUESTIONS:
!  -is the periodic correlation function correct in the program
!  -related: is it okay to take only the log of absolute values in corrln
!  Salman says we need to normalize the correlation function.  Why isn't
!  this being done already?  Although it is clear that the corrln will never
!  operate correctly if we're getting negative values for corr., since i take
!  the abs value before taking the log

!------------------------------------------------------------------------------
! lame6.hpf
! Matt Holmes, LANL T-8
! Based on the program qes.hpf written by Salman Habib.
! This program evolves a field on a lattice using the lame potential.  
! The lame potential is periodic, formed from the square of an elliptical 
! function. Its special properties of interest to us is that the equivalent 
! quantum mechanical schrodinger equation (obtained via the transfer integral 
! method) is exactly solvable. This will allow us to compare the numerical 
! correlation length and pdf calculated by this program with the analytic 
! results.

! Log of modifications from the previous qes7.hpf: 
! A) Include subroutine to calculate sn(x,k),cn(x,k),dn(x,k).  Subroutine 
!   comes from Thompson.  
! B) Modify subroutine "pot" to print out the lame potential (uses sn).
! C) Modify subroutine "rk2" to implement the lame force (uses sn,dn,cn).
! D) Modify subroutine "pot" to print out force also (if desired).
! E) Calculate periodic correlations and correlation length
! F) Take log only of absolute values in corrln
! G) Subroutine perpdf will periodify the pdfsym, output pdf_sum
! H) Calculate periodic correlation using a cosine correlator
! I) Increase the number of output files in order to "restart" the simulation.

! -----------------------------------------------------------------------------
! LIST OF OUTPUT FILES ("!" signifies that file output is commented out):
! NOTE: We care most about 85 and 96, with the pdf and the correlation length.
! fort.24 x,sn,cn,dn
! fort.25  x,potential
! fort.26  x,force
! fort.60  itime, penergy, tenergy, kenergy
! fort.64 EQUIL: xlatticepoint, averaged field-field correlation, 
!                fld^2-fld^2 correlation
! fort.65 FINAL: xlatticepoint, averaged field-field correlation, 
!                fld^2-fld^2 correlation
! fort.80 EQUIL: latticebin, pdf, pdfsym
! fort.81 FINAL: latticebin, pdf, pdfsym 
! fort.85 FINAL:  phi,perpdf (for "periodic pdf")
! fort.90  grnd state energy values
! fort.96  xlatticepoint, scaled log field-field correlation, same for squared corr.
! fort.195 every data set of the field-field correlations, at each checkpoint
! fort.98 x,v
! fort.99 x,phi
! fort.311 phi
! fort.312 v 
! fort.313 sumc1
! fort.314 sumc2
! fort.315 prevtime
! fort.316 pdf

       PROGRAM lame

!-----Langevin solver for double well QES field theory.
!-----Salman Habib, April 27, 1997. (uses Grant's stochastic RK2)

      implicit none

!-----DEFINE INPUT VARIABLES-------------------------------------------------

      integer :: nx,np,npdf,pdstart,integer,pdfint
      integer :: interval,tottime,newdata,stime,naves

      real :: eps,theta,etabar,length,xxl,n,k

      print *, "program is running..."
      print *, "opening input file now"

	open(33,file='lame6.in')
	
	read(33,*)newdata !A logial saying whether to continue an old sim.
	read(33,*)theta   !Temperature
	read(33,*)etabar  !Viscosity coefficient

	read(33,*)nx      !Number of spatial lattice points
	read(33,*)np      !Number of correlation function lattice points.
	read(33,*)npdf    !Number of histrogramming bins for the dist function

	read(33,*)eps     !Delta t
	read(33,*)length  !Delta x

	read(33,*)pdfint  !Normalizing factor for the dist function
	read(33,*)pdstart !Time at which to begin taking data (after equil)
	read(33,*)interval!Amount of time to evolve between diagnostic checks
	read(33,*)stime   !Start time
	read(33,*)tottime !Total simulation time
	read(33,*)xxl     !Used to determine range of potential to plot
	read(33,*)n	  !A parameter in the lame potential
	read(33,*)k	  !The "modulus" for the elliptical function in the 
			  ! lame potential
	read(33,*)naves   !The number of images of the shifted pdf to average
	close(33)

	print *, "Read input file."

!-----CALL SOLVE (Master Program)--------------------------------------------

      call solve(newdata,theta,etabar,nx,np,npdf,eps,length,  
     &     pdfint,pdstart,interval,stime,tottime,n,k,xxl,naves)

      end PROGRAM lame

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!-----SUBROUTINE  SOLVE (Master Program)-------------------------------------

      SUBROUTINE solve(newdata,theta,etabar,nx,np,npdf,eps,length, 
     &      pdfint,pdstart,interval,stime,tottime,n,k,xxl,naves)

      implicit none

!-----FRONT END VARIABLES----------------------------------------------------
! Dummy variables passed to subroutine:
      integer,intent(in) :: newdata,nx,np,npdf,pdfint,pdstart,
     &   interval,stime,tottime,naves
      real,intent(in) :: theta,etabar,eps,length,n,k,xxl

! Local variables:
      integer :: i,itime,otime,ftime,prevtime 
      real :: int,spacing

      integer,dimension(npdf) :: pdf,add
!hpf$ distribute pdf(block)
!hpf$ align (:) with pdf(:) :: add

      real,dimension(npdf) :: pdfsym, pdf_sum
!hpf$ distribute pdfsym(block)
!hpf$ align (:) with pdfsym(:) :: pdf_sum

      real,dimension(nx) :: phi,noise,v
!hpf$ distribute phi(block)
!hpf$ align (:) with phi(:) :: noise,v

      real,dimension(np) :: sumc1,sumc2
      real,dimension(np) :: c1add,c2add,sumc1w,sumc2w
!hpf$ distribute sumc1(block)
!hpf$ align (:) with sumc1(:) :: sumc2,c1add,c2add,sumc1w,sumc2w

!-------F90 RANDOM NUMBER STUFF-----------------------------------------------

! why the hell is this needed?  we don't change the value of the seed, and 
! we don't use this array anywhere else
!      integer,dimension(:),allocatable :: seed
      integer seed_size
      call random_seed(SIZE=seed_size)
      print *, "seed_size is ", seed_size
!      allocate (seed(seed_size))
!      call random_seed(GET=seed(1:seed_size))

! if we want to set the seed, we should read "seed" from the input file, and
! use this next line
!	random_seed(PUT=seed)
! where seed should be an integer larger than seed_size

	print *, "Allocated main arrays, set random seed size"

!-------PRINT OUT POTENTIAL IF NEEDED-----------------------------------------
	print *, "Printing the potential to file 'fort.25'"
        call pot(n,k,xxl,spacing)

!-------READ IN INPUT FILES IF NEEDED-----------------------------------------
! The information we need to continue the simulation is phi and v. 
! The information we need to continue the diagnostics are sumc1,sumc2,pdf;
! also, we need to keep track of the number of data sets to average
! over; I stored this as "prevtime".

	if (newdata.eq.1) then
           print *, "Continuing old simulation, reading files."

		open(311,file='fort.311',status='unknown')
   		read(311,*)phi
       		close(311)

		open(312,file='fort.312',status='unknown')
   		read(312,*)v
       		close(312)

                open(313,file='fort.313',status='unknown')
                read(313,*)sumc1
                close(313)

                open(314,file='fort.314',status='unknown')
                read(314,*)sumc2
                close(314)

                open(315,file='fort.315',status='unknown')
                read(315,*)prevtime
                close(315)

                open(316,file='fort.316',status='unknown')
                read(316,*)pdf
                close(316)

	else if (newdata.eq.0) then
!-------INITIALIZE PHI IF FILES ARE NOT READ IN-------------------------------

	phi=0
	v=0
        sumc1=0.0
        sumc2=0.0
        pdf=0.0
        prevtime=0
	end if

!-------INITIALIZATIONS-------------------------------------------------------

	pdfsym=0
	pdf_sum=0
        otime=0
	c1add=0
        c2add=0

!-------PERFORM EQUILIBRATION ONLY FOR NEW SIMULATION
        if (newdata.eq.0) then

	print *, "Starting equilibration loop now"

!-------LOOPING STARTS (INITIAL LOOP)-----------------------------------------

	do itime=stime,pdstart

!-------2nd order RK Langevin step

        call random_number(noise)
	where(noise.gt.5.0/6.0)noise=sqrt(6.0*etabar*theta/length/eps)
	where(noise.le.1.0/6.0)noise=-sqrt(6.0*etabar*theta/length/eps)
	where(abs(noise).le.1.0)noise=0.0
        call rk2(phi,v,noise,length,etabar,eps,n,k,nx)

        enddo

	print *, "Finished equil. loop, performing equil diagnostics."

!-------Call diagnostics/analysis routines

        itime=pdstart

	print *, "Calling unpdf:"

        call unpdf(phi,add,nx,pdfint,npdf)
        pdf=pdf+add
!UNCOMMENT for hpf: 
        forall(i=1:npdf) pdfsym(i)=(pdf(i)+pdf(npdf+1-i))/2.0
!UNCOMMENT for f90:
!        do i=1,npdf
!           pdfsym(i)=(pdf(i)+pdf(npdf+1-i))/2.0
!        end do

        pdfsym=pdfsym/(1.0*(prevtime+1))
        int=(sum(pdfsym))/(1.0*pdfint)
        pdfsym=pdfsym/int

!        call energy(phi,v,length,n,k,nx,itime)
!        call grndenergy(pdfsym,npdf,pdfint,itime,n,k,theta)

	print *, "Post equilibration, calling corrbf1"

        call corrbf1(phi,c1add,nx,np,spacing)
        sumc1=sumc1+c1add
        sumc1w=sumc1/(1.0*(prevtime+1)) !Average the correlations.       

	print *, "Post equilibration, calling corrbf2"

        call corrbf2(phi,c2add,nx,np,spacing)
        sumc2=sumc2+c2add
        sumc2w=sumc2/(1.0*(prevtime+1)) !Average the correlations.       

        print *, "Calling perpdf"

	call perpdf(npdf,pdfsym,pdf_sum,spacing,pdfint,naves)

	print *, "Finished equilibrium loop, performing diagnostics..."

        call printi(pdf,pdfsym,npdf,pdfint,np,sumc1w,sumc2w)

        end if

	print *, "starting main loop"

!-------LOOPING STARTS (NESTED LOOP)------------------------------------------
!-------Outer Loop: Carry out diagnostics/analysis at end of each step

	do otime=1,(tottime-pdstart)/interval

	print *, "One iteration of main loop going now..."

!-------Second Loop: Time step through without interruptions

        do ftime=1,interval

        call random_number(noise)
	where(noise.gt.5.0/6.0)noise=sqrt(6.0*etabar*theta/length/eps)
	where(noise.le.1.0/6.0)noise=-sqrt(6.0*etabar*theta/length/eps)
	where(abs(noise).le.1.0)noise=0.0
	
        call rk2(phi,v,noise,length,etabar,eps,n,k,nx)

        end do

!-------Call diagnostics/analysis routines

        itime=pdstart+otime*interval

	print *, "In main loop, calling unpdf"

        call unpdf(phi,add,nx,pdfint,npdf)

        pdf=pdf+add
! UNCOMMENT for hpf: 
        forall(i=1:npdf) pdfsym(i)=(pdf(i)+pdf(npdf+1-i))/2.0
! UNCOMMENT for f90:
!        do i=1,npdf
!            pdfsym(i)=(pdf(i)+pdf(npdf+1-i))/2.0
!        end do

        pdfsym=pdfsym/(1.0*(otime+prevtime+1)) !We have done otime  diag. in 
                            !main loop, 1 diag. in equil. loop, and prevtime 
                            !diag. in prev simulation= otime+prevtime+1
        int=(sum(pdfsym))/(1.0*pdfint)
        pdfsym=pdfsym/int

!        call energy(phi,v,length,n,k,nx,itime)
!        call grndenergy(pdfsym,npdf,pdfint,itime,n,k,theta)

	print *, "in main loop, calling corrbf1"

        call corrbf1(phi,c1add,nx,np,spacing)
        sumc1=sumc1+c1add

	print *, "in main loop, calling corrbf2"

        call corrbf2(phi,c2add,nx,np,spacing)
        sumc2=sumc2+c2add

	end do

	print *, "Finished main loop, performing final diagnostics..."

!-------Final computations/outputs

        sumc1w=sumc1/(1.0*(otime+prevtime+1)) !Average the correlations.
        sumc2w=sumc2/(1.0*(otime+prevtime+1))

        prevtime=prevtime+otime   !Record the number of data sets we have 
                                  !averaged over,to be used in continuing 
                                  !the simulation.
	print *, "otime= ",otime
	print *, "prevtime= ",prevtime

        print *, "Calling perpdf..."
	call perpdf(npdf,pdfsym,pdf_sum,spacing,pdfint,naves)

	print *, "Calling corrln (fort.96)..."
        call corrln(sumc1w,sumc2w,length,np)

	print *, "Calling printj (fort.65, fort.81,fort.85)..."
        call printj(pdf,pdfsym,pdf_sum,npdf,pdfint,np,sumc1w,sumc2w)

	print *, "Calling printfinal (fort.98,99,311,312,313,314,315,316)..."
        call printfinal(phi,v,nx,np,npdf,sumc1,sumc2,pdf,prevtime)

        end SUBROUTINE solve

!-------SUBROUTINE PRINTI-----------------------------------------------------

        SUBROUTINE printi(pdf,pdfsym,npdf,pdfint,np,sumc1w,  
     &     sumc2w)

        implicit none

! Dummy variables:
        integer, intent(in) ::  pdfint,npdf,np

        integer, dimension(npdf) :: pdf
!hpf$ distribute pdf(block)

        real, dimension(npdf) :: pdfsym,xxx        
!hpf$ distribute pdfsym(block)
!hpf align (:) with pdfsym(:) :: xxx

        real, dimension(np) :: sumc1w,sumc2w
!hpf$ distribute sumc1w(block)
!hpf align (:) with sumc1w(:) :: sumc2w

! Local variables:
        integer :: i, nxm
        real :: xxmin

        xxmin=-0.50*(1.0*npdf-1.0)/(1.0*pdfint)
!UNCOMMENT for hpf: 
        forall(i=1:npdf) xxx(i)=xxmin + (i-1)/(1.0*pdfint)
!UNCOMMENT for f90:
!        do i=1,npdf
!            xxx(i)=xxmin + (i-1)/(1.0*pdfint)
!        end do

        do i=1,np
        nxm=i-1
        write(64,*)nxm,sumc1w(i),sumc2w(i)
        enddo

        do i=1,npdf
        write(80,*)xxx(i),pdf(i),pdfsym(i)
        enddo

	end SUBROUTINE printi

!-------SUBROUTINE PRINTJ-----------------------------------------------------

        SUBROUTINE printj(pdf,pdfsym,pdf_sum,npdf,pdfint,np,sumc1w,  
     &     sumc2w)

        implicit none

! Dummy variables:
        integer, intent(in) ::  pdfint,npdf,np

        integer, dimension(npdf) :: pdf
!hpf$ distribute pdf(block)

        real, dimension(npdf) :: pdfsym,pdf_sum,xxx        
!hpf$ distribute pdfsym(block)
!hpf$ align (:) with pdfsym(:) :: xxx,pdf_sum

        real, dimension(np) :: sumc1w,sumc2w
!hpf$ distribute sumc1w(block)
!hpf$ align (:) with sumc1w(:) :: sumc2w

! Local variables:
        integer :: i, nxm
        real :: xxmin

        xxmin=-0.50*(1.0*npdf-1.0)/(1.0*pdfint)
!UNCOMMENT for hpf: 
        forall(i=1:npdf) xxx(i)=xxmin + (i-1)/(1.0*pdfint)
!UNCOMMENT for f90:
!        do i=1,npdf
!            xxx(i)=xxmin + (i-1)/(1.0*pdfint)
!        end do

        do i=1,np
        nxm=i-1
        write(65,*)nxm,sumc1w(i),sumc2w(i)
        enddo

        do i=1,npdf
        write(81,*)xxx(i),pdf(i),pdfsym(i)
        enddo

	do i=1,npdf
	write(85,*)xxx(i),pdf_sum(i)
	enddo

	end SUBROUTINE printj

!-------SUBROUTINE PRINTFINAL--------------------------------------------------

        SUBROUTINE printfinal(phi,v,nx,np,npdf,sumc1,sumc2,pdf,prevtime)

        implicit none
! Dummy variables:
        integer, intent(in):: nx,np,npdf
        real, dimension(nx), intent(in) :: phi,v
        real,dimension(np),intent(in) :: sumc1,sumc2
        integer,dimension(npdf),intent(in) :: pdf
        integer,intent(in) :: prevtime

! Local variables:
        integer :: i

!hpf$ distribute phi(block)
!hpf$ align (:) with phi(:) :: v

!hpf$ distribute sumc1(block)
!hpf$ align (:) with sumc1(:) :: sumc2

!hpf$ distribute pdf(block)

! This format is for plotting:
        do i=1,nx
        write(99,*)i,phi(i)
        write(98,*)i,v(i)
        enddo

! The following output files ("300"series) are for continuing another sim.
        write(311,*)phi
        write(312,*)v
        write(313,*)sumc1
        write(314,*)sumc2
        write(315,*)prevtime
        write(316,*)pdf

	end SUBROUTINE printfinal

!-------SUBROUTINE GRNDENERGY-------------------------------------------------

        SUBROUTINE grndenergy(pdfsym,npdf,pdfint,itime,n,k,theta)

        implicit none

! Dummy variables:
        integer,intent(in) ::  pdfint,npdf,itime
        real,intent(in) :: n,k,theta
        real,dimension(npdf),intent(in) :: pdfsym

! Local variables:
        integer :: i
        real :: int,xxmin,kep,vep,tep
        real,dimension(npdf) :: xxx,psip
	real,dimension(npdf) :: sn,cn,dn  !Holds values for elliptic functions

!hpf$ distribute pdfsym(block)
!hpf$ align (:) with pdfsym(:) :: xxx,psip,sn,cn,dn

        xxmin=-0.50*(1.0*npdf-1.0)/(1.0*pdfint)
!UNCOMMENT for hpf: 
        forall(i=1:npdf) xxx(i)=xxmin + (i-1)/(1.0*pdfint)
!UNCOMMENT for f90:
!        do i=1,npdf
!            xxx(i)=xxmin + (i-1)/(1.0*pdfint)
!        end do
        psip=sqrt(pdfsym)
        kep=-sum(psip*(eoshift(psip,-1,0.0,1)+eoshift(psip,1,0.0,1)-2.0*  
     &        psip)) *theta*theta*pdfint*0.50
	call sncndn(xxx,k,npdf,sn,cn,dn)
        vep=(sum(pdfsym*(n*(n+1)*k*k*sn*sn)))/(1.0*pdfint)
        tep=kep+vep
        write(90,*)itime,tep,vep,kep

        end SUBROUTINE grndenergy

!-------SUBROUTINE CORRLN-------------------------------------------------

        SUBROUTINE corrln(sumc1w,sumc2w,length,np)

        implicit none

! Dummy variables:
        integer,intent(in) :: np
        real, intent(in) :: length
        real,dimension(np),intent(in) :: sumc1w,sumc2w
        
! Local variables:
        integer :: i
        real,dimension(np) :: logc1,logc2
        real,dimension(np-1) :: gamma1,gamma2

!hpf$ distribute sumc1w(block)
!hpf$ align (:) with sumc1w(:) :: sumc2w,logc1,logc2

!hpf$ distribute gamma1(block)
!hpf$ align (:) with gamma1(:) :: gamma2

        logc1=log(abs(sumc1w))
        logc2=log(abs(sumc2w))

!UNCOMMENT for hpf: 
        forall(i=1:np-1) gamma1(i)=(logc1(i+1)-logc1(i))/length
!UNCOMMENT for hpf: 
        forall(i=1:np-1) gamma2(i)=(logc2(i+1)-logc2(i))/length
!UNCOMMENT for f90:
!        do i=1,np-1
!            gamma1(i)=(logc1(i+1)-logc1(i))/length
!        end do
!UNCOMMENT for f90:
!        do i=1,np-1
!            gamma2(i)=(logc2(i+1)-logc2(i))/length
!        end do

        do i=1,np-1
        write(96,*)i,gamma1(i),gamma2(i)
        enddo

        end SUBROUTINE corrln

!-------SUBROUTINE POT--------------------------------------------------------
! Our potential is the lame potential. It is periodic and defined with the
! elliptical function sn.  The parameters for this potential are n and k.
! The force is the derivative of the potential, and uses the functions sn, dn,
! and cn.
! This function also calculates the period of the potential for use in the 
! per_pdf subroutine.
 
        SUBROUTINE pot(n,k,xxl,spacing)

        implicit none

! Dummy variables:
        real,intent(in) :: n,k,xxl     !xxl gives the range of the pot to print
	real, intent(inout) :: spacing !spacing is the period of the potential

! Local variables:
        integer :: i,count
	real :: current,value
        real :: dx,xmin
        real,dimension(1024) :: potential,xx
	real,dimension(1024) :: sn,cn,dn
	real,dimension(1024) :: force

!hpf$ distribute potential(block)
!hpf$ align (:) with potential(:) :: xx,sn,cn,dn
!hpf$ align (:) with potential(:) :: force
  
        xmin=-0.50*xxl
        dx=xxl/(1.0*1024)
!UNCOMMENT for hpf: 
        forall(i=1:1024) xx(i)=xmin+(i-1)*dx
	call sncndn(xx,k,1024,sn,cn,dn)

!UNCOMMENT for hpf: 
        potential=n*(n+1)*k*k*sn*sn
	force=-2.0*n*(n+1)*k*k*sn*cn*dn

!UNCOMMENT for f90:
!        do i=1,1024
!            xx(i)=xmin+(i-1)*dx
!        end do
!UNCOMMENT for f90:
!        do i=1,1024
!            potential(i)=n*(n+1)*k*k*sn(i)*sn(i)
!	     force(i)=-2.0*n*(n+1)*k*k*sn(i)*cn(i)*dn(i)
!        end do

! Here we calculate the spacing between minima of the potential, to be used
! to periodify the pdf.  We do it with the force, by starting at 0 and finding
! the next zero of the force.  Then we double it to obtain the period of the pot.

        count=512       !Start at the middle of the array, where xx(512)=0
	current=0
	value=0
	spacing=0

   	do 
	count=count+1	!will be used as an array indice (start at second
                        !location in the array)
	current=force(count)	!Get next value in the array
	if(current > 0) exit
	end do

!checks:
	print *, "count is", count
	print *, "the first value greater than zero is", current
	print *, "the just previous value was", force(count-1)
	print *, "we are averaging", xx(count-1), " and", xx(count)

	value=(xx(count-1)+xx(count))/2 !use the average of the location just 
					!before the zero and just after.
	spacing=2*value         !The location of the first zero crossing is only
				!half the period of the potential.
	print *, "value is", value
	print *, "spacing is", spacing

        do i=1,1024
        write(25,*)xx(i),potential(i)
	write(26,*)xx(i),force(i)
	write(24,*)xx(i),sn(i),cn(i),dn(i)  !Write out the elliptic functions
        enddo

        end SUBROUTINE pot

!-------SUBROUTINE CORRBF1----------------------------------------------------

        SUBROUTINE corrbf1(phi,c1add,nx,np,spacing)
! Calculates a periodic field-field correlation function:
! C(u)=<|cos((pi/2)(1/min_location)(u(x)-u(0)))|>
! For this to be appropriately periodic, the values of the field must be 
!   scaled by pi/2; this is being achieved within the correlation calculation
!   itself.

        implicit none

! Dummy variables:
        integer,intent(in) :: nx,np
	real,intent(in) :: spacing
        real,dimension(np),intent(inout) :: c1add
        real,dimension(nx),intent(in) :: phi
!hpf$ distribute phi(block)

! Local variables:
        real,parameter :: pi=3.1415926, piby2=pi/2.0
        integer :: nxm,i  !nxm is used as a dummy index when writing to the file, see below

        real :: phiavg,phi2avg
        real,dimension(np) :: cphi

!hpf$ distribute c1add(block)
!hpf$ align (:) with c1add(:) :: cphi

        cphi=0

        phiavg=sum(phi)/(1.0*nx)   !averages will be used for normalization
        phi2avg=sum(phi*phi)/(1.0*nx)

        do i=1,np
! Calculate the correlation
        cphi(i)=sum(abs(cos((piby2/spacing)*(cshift(phi,i-1)-phi))))
     $          /(1.0*nx)
        enddo

! NOT SURE WHAT NORMALIZATION TO USE NOW?
!        c1add=(cphi-phiavg*phiavg)/(phi2avg-phiavg*phiavg)  !a normalization 

        c1add=cphi

        do i=1,np
        nxm=i-1
        write(195,*)nxm,c1add(i)
        enddo

        end SUBROUTINE corrbf1

!-------SUBROUTINE CORRBF2----------------------------------------------------

        SUBROUTINE corrbf2(phi,c2add,nx,np,spacing)

        implicit none

! Dummy variables:
        integer,intent(in) :: nx,np
	real,intent(in) :: spacing
        real,dimension(np),intent(inout) :: c2add
        real,dimension(nx),intent(in) :: phi
!hpf$ distribute phi(block)

! Local variables:
        real,parameter :: pi=3.1415926, piby2=pi/2.0
        integer :: nxm,i
        real :: phi2avg
        real,dimension(np) :: cphi2r,cphi2
!hpf$ distribute c2add(block)
!hpf$ align (:) with c2add(:) :: cphi2,cphi2r

        cphi2=0.0
        phi2avg=sum(phi*phi)/(1.0*nx)

! INCOMPLETE: how do we alter this to make it appropriately periodic?
        do i=1,np
        cphi2(i)=sum(abs(cos((piby2/spacing)*(cshift(phi,i-1)*
     $          cshift(phi,i-1)-phi*phi))))/(1.0*nx)
        enddo

!A NORMALIZATION I DON'T UNDERSTAND:
!        cphi2r=cphi2-phi2avg**2
         c2add=cphi2r-phi2avg*phi2avg 

!        do i=1,np
!        nxm=i-1
!        write(196,*)nxm,c2add(i)
!        enddo

        end SUBROUTINE corrbf2
       
!-------SUBROUTINE UNPDF------------------------------------------------------
 
        SUBROUTINE unpdf(phi,add,nx,pdfint,npdf)

        implicit none

! Dummy variables:
        integer,intent(in) ::  nx,pdfint,npdf
        integer,dimension(npdf),intent(inout) :: add
!hpf$ distribute add(block)
        real,dimension(nx) :: phi

! Local variables:
        integer :: i,nco
        logical,dimension(nx) :: mmsskk
!hpf$ distribute mmsskk(block)
        real,dimension(nx) :: iphi
!hpf$ distribute phi(block)
!hpf$ align (:) with phi(:) :: iphi

        iphi=1.0*pdfint*phi

        do i=-(npdf-1)/2,(npdf-1)/2

        mmsskk=.false.
        where(i.le.iphi.and.iphi.lt.i+1)mmsskk=.true.
        nco=count(mmsskk)
        add((npdf-1)/2+1+i)=nco

	enddo

        end SUBROUTINE unpdf

!-------SUBROUTINE ENERGY-----------------------------------------------------

        SUBROUTINE energy(phi,v,length,n,k,nx,itime)

        implicit none

! Dummy variables:
        integer,intent(in) ::  nx,itime
        real,intent(in) ::  length,n,k
        real,dimension(nx),intent(in) :: phi
!hpf$ distribute phi(block)

! Local variables:
        real :: kenergy,penergy,tenergy
        real,dimension(nx) :: gr,v
	real,dimension(nx) :: sn,cn,dn
!hpf$ align (:) with phi(:) :: gr,v,sn,cn,dn

        gr=(cshift(phi,1,1)-phi)/length
	call sncndn(phi,k,nx,sn,cn,dn)
        penergy=length*sum(.50*gr*gr+n*(n+1)*k*k*sn*sn)
        kenergy=.50*length*sum(v*v)
        tenergy=kenergy+penergy
  
        penergy=penergy/(nx*length)
        kenergy=kenergy/(nx*length)
        tenergy=tenergy/(nx*length)

        write(60,*)itime,tenergy,penergy,kenergy

        end SUBROUTINE energy

!-------SUBROUTINE RK2--------------------------------------------------------

        SUBROUTINE rk2(phi,v,noise,length,etabar,eps,n,k,nx)
        implicit none

! Dummy variables:
        integer,intent(in) :: nx
        real,intent(in) :: length,etabar,eps,n,k
        real,dimension(nx),intent(in) :: noise
        real,dimension(nx),intent(inout)::phi,v

! Local variables:
	real,dimension(nx) :: sn,cn,dn !arrays to hold the values of elliptic 
				!functions for each value of the phi array
        real,dimension(nx) ::pibefore,lapphi, phibefore,tv1,tv2,tphi1
!hpf$ align (:) with phi(:) :: sn,cn,dn
!hpf$ align (:) with phi(:) :: pibefore,phibefore,lapphi,v,noise,tv1,tv2,tphi1

        pibefore=v
        phibefore=phi

        lapphi=(cshift(phi,-1)+cshift(phi,1)-2.0*phi)/length/length
! Derivative of the potential is:
! dV/du = 2 n(n+1)k^2 cn(u,k)dn(u,k)sn(u,k)
	call sncndn(phi,k,nx,sn,cn,dn)  !Calculate elliptic functions for all
					! values of phi
        tv1=lapphi - etabar*v + noise     
     &             - 2.0*n*(n+1)*k*k*sn*cn*dn
        tphi1=v
        phi=phibefore+eps*tphi1
        v=v+eps*tv1

        lapphi=(cshift(phi,-1)+cshift(phi,1)-2.0*phi)/length/length
! changed potential here also...
	call sncndn(phi,k,nx,sn,cn,dn)  !Calculate elliptic functions
        tv2=lapphi - etabar*v + noise   
     &             - 2.0*n*(n+1)*k*k*sn*cn*dn
        v=pibefore+eps*0.5*(tv1+tv2)
        phi=phibefore+eps*0.5*(tphi1+v)

        end SUBROUTINE rk2

!-------SUBROUTINE sncndn------------------------------------------------------

! sn,cn, and dn are jacobian elliptic functions
! In our case, these functions are used to calculate the potential and force. 
! For us, "u" will be the value of the field, phi, while "k" determines the 
! shape of the potential.  k=0 reduces to a sin curve, while k=1 reduces to 
! the tanh curve. 
	
	SUBROUTINE sncndn(u,k,nx,sn,cn,dn)
! Jacobi elliptic functions by AGM method
        implicit none

! Dummy variables 
	integer, intent(in):: nx
	real, dimension(nx), intent(in):: u
	real, dimension(nx), intent(out):: sn,cn,dn   !Contains the results.
	real, intent(in):: k 

!hpf$ align (:) with u(:) :: sn,cn,dn 

! Local variables
	integer :: i     !Loop variable
	real :: eps=10E-8   !Determines the precision 
	real, dimension(0:9):: a,b,c
! we need a 9 element array for each element of the array u, in order to 
!   calculate an intermediate value.  We create a 2 dimensional array.
	real, dimension(0:9,nx):: eeks

!hpf& align eeks(*,:) with u(:)

	integer :: n,max=10

	a(0)=1.0
	b(0)=sqrt(1.0-k*k)
	c(0)=sqrt(1.0*k*k)

	n=0

	do 
	  if(.NOT.(abs(c(n)) > eps)) exit
	  a(n+1) = 0.5*(a(n)+b(n))
	  b(n+1) = sqrt(a(n)*b(n))
	  c(n+1) = 0.5*(a(n)-b(n))

	  if(n >= max-2) then
	    write(*,*) "!!In JacobiSn n=", n, " too large"
	    write(*,*) "Increase max or decrease eps"
	    c(max-1) = 0.0
	  endif

	  n=n+1
	end do

	n=n-1

	forall(i=1:nx)eeks(n,i)=(2**n)*a(n)*u(i)

	do while (n>0)
	  forall(i=1:nx)eeks(n-1,i) = 0.5*(eeks(n,i)+asin(c(n)*
     &      sin(eeks(n,i))/a(n)))
	  n=n-1
	end do

	sn(:) = sin(eeks(0,:))
	cn(:) = cos(eeks(0,:))
	dn(:) = cn/cos(eeks(1,:)-eeks(0,:))
	
	end SUBROUTINE sncndn

       
!-------SUBROUTINE PERPDF------------------------------------------------------
! The purpose of this subroutine is to "periodify" the numerical pdf.

        SUBROUTINE perpdf(npdf,pdfsym,pdf_sum,spacing,pdfint,naves)

        implicit none

! Dummy variables:
! "naves" is number of averages
        real,intent(in) :: spacing
        integer,intent(in) ::  naves,npdf,pdfint
        real,dimension(npdf),intent(inout) :: pdf_sum
        real,dimension(npdf), intent(in) :: pdfsym        
!hpf$ distribute pdfsym(block)
!hpf$ align (:) with pdfsym(:) :: pdf_sum

! Local variables:
        integer :: k,lspacing

        lspacing=0

        print *, "spacing is", spacing
        print *, "pdfint is", pdfint

! The scaling between phi and the npdf sized lattice is given by pdfint
	lspacing=int(spacing*pdfint) !cast it to an integer

! check
	print *, "lspacing is ", lspacing

!question: do I have to perform an explicit summing over the indices of the 
!array?

        print *, "performing the periodification process promptly"

	do k=1,naves
	  pdf_sum=pdf_sum+
     $          eoshift(pdfsym,((-((naves-1)/2)-1)+k)*lspacing)
	enddo

! Above, we started the averaging with the left-most shifted pdf, then 
! worked our way right until adding all of the pdfs

! Perform the average over the number of images summed up
	pdf_sum=pdf_sum/naves

        end SUBROUTINE perpdf