title: debugging old fortran code
body: <p>I have an old Fortran code for calculation of Lyapunov exponent which I tried converting to modern Fortran syntax.</p>
  
  <pre><code>PROGRAM ODE 
  integer, PARAMETER :: N=3 
  integer, PARAMETER :: NN=12 
  EXTERNAL FCN 
  DIMENSION Y(NN),ZNORM(N),GSC(N),CUM(N),C(24),W(NN,9) 
  Y(1) = 10.0 
  Y(2) = 1.0 
  Y(3) = 0.0 
  ! INITIAL CONDITIONS FOR LINEAR SYSTEM (ORTHONORMAL FRAME) 
  DO 10 I = N+1,NN 
  Y(I) = 0.0 
  10 CONTINUE 
  DO 20 I = 1,N 
  Y((N+1)*I) = 1.0 
  CUM(I) = 0.0 
  20 CONTINUE 
  ! INTEGRATION TOLERANCE, # OF INTEGRATION STEPS, 
  ! TIME PER STEP, AND I/O RATE 
  write (*,*) "TOL, NSTEP, STPSZE, IO ?" 
  read  (*,*) TOL, NSTEP, STPSZE, IO 
  ! INITIALIZATION FOR INTEGRATOR 
  NEQ = NN 
  X=0.0
  IND = 1
  DO 100 I = 1,NSTEP 
  XEND = STPSZE*FLOAT(I) 
  ! CALL ANY ODE INTEGRATOR - THIS IS AN LMSL ROUTINE 
  CALL DVERK (NEQ,FCN,X,Y,XEND,TOL, IND,C,NEQ,W,IER) 
  ! CONSTRUCT A NEW ORTHONORMAL BASIS BY GRAM-SCHMIDT METHOD 
  ! NORMALIZE FIRST VECTOR 
  ZNORM(1) = 0.0 
  DO 30 J = 1,N 
  ZNORM(1) = ZNORM(1)+Y(N*J+1)**2 
  30 CONTINUE 
  ZNORM(1) = SQRT(ZNORM(1)) 
  DO 40 J = 1,N 
  Y(N*J+1) = Y(N*J+1)/ZNORM(1) 
  40 CONTINUE 
  ! GENERATE THE NEW ORTHONORMAL SET OF VECTORS. 
  DO 80 J = 2,N 
  ! GENERATE J-1 GSR COEFFICIENTS. 
  DO 50 K = l,(J-l) 
  GSC(K) = 0.0 
  DO 50 L = 1,N 
  GSC(K) = GSC(K)+Y(N*L+J)*Y(N*L+K) 
  50 CONTINUE 
  ! CONSTRUCT A NEW VECTOR. 
  DO 60 K = 1,N 
  DO 60 L = l,(J-l) 
  Y(N*K+J) = Y(N*K+J)-GSC(L)*Y(N*K+L) 
  60 CONTINUE 
  ! CALCULATE THE VECTOR'S NORM 
  ZNORM(J) = 0.0 
  DO 70 K = I,N 
  ZNORM(J) = ZNORM(J)+Y(N*K+J)**2 
  70 CONTINUE 
  ZNORM(J) = SQRT(ZNORM(J)) 
  ! NORMALIZE THE NEW VECTOR. 
  DO 80 K = 1,N 
  Y(N*K+J) = Y(N*K+J)/ZNORM(J) 
  80 CONTINUE 
  ! UPDATE RUNNING VECTORMAGNITUDES 
  DO 90 K = 1,N 
  CUM(K) = CUM(K)+ALOG(ZNORM(K) )/ALOG(2. ) 
  90 CONTINUE 
  ! NORMALIZE EXPONENT AND PRINT EVERY IO ITERATIONS 
  IF (MOD(I,IO).EQ.0) write (*,*) X,(CUM(K)/X,K = I,N) 
  100 CONTINUE 
  CALL EXIT 
  END 
  SUBROUTINE FCN (N,X,Y,YPRIME) 
  ! USER DEFINED ROUTINE CALLED BY IMSL INTEGRATOR. 
  DIMENSION Y(12),YPRIME(12) 
  ! LORENZ EQUATIONS OF MOTION 
  YPRIME(1) = 16.*(Y(2)-Y(1)) 
  YPRIME(2) = -Y(1)*Y(3)+45.92*Y(1)-Y(2) 
  YPRIME(3) = Y(1)*Y(2)-4.*Y(3) 
  ! 3 COPIES OF LINEARIZED EQUATIONS OF MOTION. 
  DO 10 I = 0,2 
  YPRIME(4+I) = 16.*(Y(7+I)-Y(4+I)) 
  YPRIME(7+I) = (45.92-Y(3))*Y(4+I)-Y(7+I)-Y(1)*Y(10+I) 
  YPRIME(10+I) = Y(2)*Y(4+I)+Y(1)*Y(7+I)-4.*Y(10+I) 
  10 CONTINUE 
  RETURN 
  END 
  </code></pre>
  
  <p>I have debugged most of this, but I am still left with a few errors that I am unable to get around. The error log says:    </p>
  
  <pre><code>main.f95:44.14:
  
  DO 50 L = 1,N 
                1
  Warning: Obsolescent feature: Shared DO termination label 50 at (1)
  main.f95:49.18:
  
  DO 60 L = l,(J-l) 
                    1
  Warning: Obsolescent feature: Shared DO termination label 60 at (1)
  main.f95:59.14:
  
  DO 80 K = 1,N 
                1
  Warning: Obsolescent feature: Shared DO termination label 80 at (1)
  /tmp/ccfI69Sj.o: In function `MAIN__':
  main.f95:(.text+0x296): undefined reference to `dverk_'
  main.f95:(.text+0x844): undefined reference to `exit_'
  collect2: error: ld returned 1 exit status
  </code></pre>
  
  <p>Could someone please help me out in resolving the errors?
  Thanks.</p>