Clicky

Fortran Wiki
Fixed form layout

Below is an example of a program written in fixed form layout, inherited from punched cards, where column one contains *, !, or C for comments, column six may contain a continuation character, and executable statements begin at column seven. Fortran 90 introduced the Free form layout.

     program sphere

c ********** Distribution of eigenvalues in a sphere of radius R
c ********** n atoms, all distances in units of 1/k

      complex*16 c1
      real*8 pi
      integer n, nreal, lambda
      integer lwork, ldvl, ldvr
      parameter(c1 = (0.0,1.0))
      parameter(n = 1000)
      parameter(lwork = 4*n)
      parameter(ldvl = 1)
      parameter(ldvr = 1)
      parameter(pi = 3.141592653589793238462643383279502884197)
      parameter(nreal = 2000)
      parameter(nlambda = 300)

      real*8 akR, akR2, rho, beta
      real*8 tmp1, tmp2, tmp3
      real*8 x1(1:n), x2(1:n), x3(1:n)
      real*8 p2d(1:nlambda,1:nlambda)
	real*8 pre(1:nlambda), pim(1:nlambda)
      real*8 lambdare1, lambdare2, dlambdare
      real*8 lambdaim1, lambdaim2, dlambdaim
	real*8 rwork(1:2*n)

      complex*16 a(1:n,1:n)
      complex*16 vl1(1:ldvl, 1:n), vr1(1:ldvr, 1:n)
      complex*16 vl(1:n, 1:n), vr(1:n, 1:n)
      complex*16 w(1:n)
      complex*16 work(1:lwork)
	complex*16 ctmp1, ctmp2, ctmp3, ctmp4

      integer info
      integer i, j, k, m, ireal, nrot

      external RandomNumber
      real*8 RandomNumber

c ********** Initialization of parameters

      rho = 10.0D+00

      akR = (2.0D+00)*pi*((3.0D+00*n/
    *       ((4.0D+00)*pi*rho))**((1.0D+00)/(3.0D+00)))
	akR2 = akR**2

c	beta = (2.8D+00)*n/((akL**2))

	 write(6,*) "N = ", n
	 write(6,*) "rho = ", rho
	 write(6,*) "kR = ", akR

       do i = 1, nlambda
         pre(i) = 0.0D+00
         pim(i) = 0.0D+00
	   do j = 1, nlambda
	     p2d(i,j) = 0.0D+00
	   end do
       end do

	lambdare1 = -5.0D+00
	lambdare2 = 5.0D+00
	dlambdare = (lambdare2 - lambdare1)/nlambda

      lambdaim1 = -1.0D+00
	lambdaim2 = 7.0D+00
	dlambdaim = (lambdaim2 - lambdaim1)/nlambda

      evmin1 = 1000.0D+00
	evmin2 = 1000.0D+00
	evmax = -1000.0D+00

	do ireal = 1, nreal

        write(6,*) "Running", ireal, " of", nreal

c ********** Initialization of the matrix gamma
        do i = 1, n
          do j = 1, n
            a(i,j) = 0.0D+00
          end do
        end do

	  do i = 1, n
	    tmp1 = (2.0D+00)*akR2
	    do while(tmp1.GT.akR2)
	      x1(i) = akR*(-1.0D+00 + (2.0D+00)*RandomNumber(idummy))
	      x2(i) = akR*(-1.0D+00 + (2.0D+00)*RandomNumber(idummy))
	      x3(i) = akR*(-1.0D+00 + (2.0D+00)*RandomNumber(idummy))
	      tmp1 = x1(i)**2 + x2(i)**2 + x3(i)**2
	    end do
	  end do

	  do i = 1, n
	    do j = 1, n
	      a(i,j) = (x1(i) - x1(j))**2 + (x2(i) - x2(j))**2 +
    *                (x3(i) - x3(j))**2
	    end do
	  end do

        do i = 1, n
          do j = 1, n
            if(i.EQ.j) then
		 a(i,j) = 0.0D+00
            else
		 tmp1 = cdsqrt(a(i,j))
		 a(i,j) = cdexp(c1*tmp1)/tmp1
c		 a(i,j) = c1*dsin(tmp1)/tmp1
	      end if
	    end do
	  end do

c ********* Using LAPACK (compile with -llapack)

c ********* Compute only the eigenvalues
c         call zgeev('N', 'N', n, a, n, w, vl1, ldvl,
c     * vr1, ldvr, work, lwork, rwork, info)

c ********* Compute both the eigenvalues and left and right eigenvectors
        call zgeev('V', 'V', n, a, n, w, vl, n,
    * vr, n, work, lwork, rwork, info)


	  if(info.NE.0) then
          write(6,*) "Error in ZGEEV: info = ", info
	    pause
	  end if

c ********* Using IMSL
c	    CALL DEVLCG (n, a, n, w)

c ********* Save all eigenvalues

         if(ireal.EQ.1) then
	     open(UNIT=12, FILE="ev.dat",
    *STATUS="UNKNOWN")
	   else
           open(UNIT=12, FILE="ev.dat",
    *STATUS="UNKNOWN", ACCESS="APPEND")
	   end if

         do i = 1, n
           write(12,*) dreal(w(i)), " ", dimag(w(i))
         end do
         close(12)

c ********* Compute and save IPR

         if(ireal.EQ.1) then
	     open(UNIT=12, FILE="ipr.dat",
    *STATUS="UNKNOWN")
	   else
           open(UNIT=12, FILE="ipr.dat",
    *STATUS="UNKNOWN", ACCESS="APPEND")
	   end if

         do i = 1, n
	     tmp1 = 0.0D+00
	     tmp2 = 0.0D+00
	     do j = 1, n
		tmp1 = tmp1 + cdabs(vr(j,i))**2
		tmp2 = tmp2 + cdabs(vr(j,i))**4
	     end do
	     tmp3 = tmp2/(tmp1*tmp1)
           write(12,*) dreal(w(i)), " ", dimag(w(i)), " ", tmp3
         end do
         close(12)

c ********* Renormalize eigenfunctions to have <L|R> = 1

	   do i = 1, n
	     ctmp1 = 0.0D+00
	     do j = 1, n
c		ctmp1 = ctmp1 + dconjg(vl(j,i))*vr(j,i)
		ctmp1 = ctmp1 + vr(j,i)**2
	     end do
	     ctmp2 = (1.0D+00)/cdsqrt(ctmp1)
	     
	     do j = 1, n
c		vl(j,i) = vl(j,i)*dconjg(ctmp2)
		vr(j,i) = vr(j,i)*ctmp2
	     end do

	   end do

c ********* Save <L|L><R|R>

         if(ireal.EQ.1) then
	     open(UNIT=12, FILE="corr.dat",
    *STATUS="UNKNOWN")
	   else
           open(UNIT=12, FILE="corr.dat",
    *STATUS="UNKNOWN", ACCESS="APPEND")
	   end if

         do i = 1, n
	     tmp1 = 0.0D+00
	     tmp2 = 0.0D+00
	     do j = 1, n
		tmp1 = tmp1 + cdabs(vr(j,i))**2
c		tmp2 = tmp2 + cdabs(vl(j,i))**2
	     end do
c	     tmp3 = tmp1*tmp2
	     tmp3 = tmp1*tmp1/n
           write(12,*) dreal(w(i)), " ", dimag(w(i)), " ", tmp3
         end do
         close(12)

c ********** Compute probability distribution

	    do i = 1, n
	       j = int((dreal(w(i)) - lambdare1)/dlambdare) + 1
            if((j.GE.1).AND.(j.LE.nlambda)) then
		     pre(j) = pre(j) + 1.0D+00
            else
              write(6,*) j, " beyond 1:", nlambda
c               pause
            end if
	    end do

	    do i = 1, n
	       j = int((dimag(w(i)) - lambdaim1)/dlambdaim) + 1
            if((j.GE.1).AND.(j.LE.nlambda)) then
		     pim(j) = pim(j) + 1.0D+00
            else
              write(6,*) j, " beyond 1:", nlambda
c               pause
            end if
	    end do

	    do i = 1, n
	      j = int((dreal(w(i)) - lambdare1)/dlambdare) + 1
		  k = int((dimag(w(i)) - lambdaim1)/dlambdaim) + 1
           if((j.GE.1).AND.(j.LE.nlambda)) then
			if((k.GE.1).AND.(k.LE.nlambda))	then
		      p2d(j,k) = p2d(j,k) + 1.0D+00
             else
			  write(6,*) k, " beyond 1:", nlambda
c			  pause
			end if
	      else
			  write(6,*) j, " beyond 1:", nlambda
c			  pause
		  end if
         end do

      end do

c ********** Normalization

	tmp1 = (1.0D+00)/(n*nreal*dlambdare)
	tmp2 = (1.0D+00)/(n*nreal*dlambdaim)
	tmp3 = (1.0D+00)/(n*nreal*dlambdare*dlambdaim)

	do i = 1, nlambda
	  pre(i) = tmp1*pre(i)
	  pim(i) = tmp2*pim(i)
	  do j = 1, nlambda
	    p2d(i,j) = tmp3*p2d(i,j)
	  end do 
     end do

c ********** Save P(lambda)

      open(UNIT=7, FILE="pre.dat", STATUS="UNKNOWN")
      do i = 1, nlambda
        tmp1 = lambdare1 + (i - 0.5D+00)*dlambdare
        write(7,*) tmp1, " ", pre(i)
      end do
      close(7)

      open(UNIT=7, FILE="pim.dat", STATUS="UNKNOWN")
      do i = 1, nlambda
        tmp1 = lambdaim1 + (i - 0.5D+00)*dlambdaim
        write(7,*) tmp1, " ", pim(i)
      end do
      close(7)

      open(UNIT=7, FILE="p2d.dat", STATUS="UNKNOWN")
      do i = 1, nlambda
        tmp1 = lambdare1 + (i - 0.5D+00)*dlambdare
	   do j = 1, nlambda
          tmp2 = lambdaim1 + (j - 0.5D+00)*dlambdaim
          write(7,*) tmp1, " ", tmp2, " ", p2d(i,j)
	   end do
      end do
      close(7)

      end

c ***************************************

      real*8 function RandomNumber(idummy)

	integer idummy
	real*8 tmp1

c	 RandomNumber = rand()

      call RANDOM_NUMBER(tmp1)
	RandomNumber = tmp1

	return
	end