Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

  • Congratulations strongm on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!

Compiler Fortran RADAU

Status
Not open for further replies.

pistola

Programmer
Oct 1, 2018
1
BR
Hi everyone, I can not compile a program here (RADAU-INSTRUCOES.for), which uses the program radau ra15, (ra27.for)
When I compile the program RADAU-INSTRUCOES.for, the following errors appear:

Code:
./RADAU-INSTRUCOES.for:229:38:

       f(1)=-aa*x(1)-bb*x(2)-0.002d0*x(3) ! f (i) are the first derivatives of x, y, z
                                      1
Warning: Array reference at (1) is out of bounds (3 > 1) in dimension 1
./RADAU-INSTRUCOES.for:229:25:

       f(1)=-aa*x(1)-bb*x(2)-0.002d0*x(3) ! f (i) are the first derivatives of x, y, z
                         1
Warning: Array reference at (1) is out of bounds (2 > 1) in dimension 1
./RADAU-INSTRUCOES.for:230:30:

       f(2)=-bb*x(2)+ff*x(3)-x(2)
                              1
Warning: Array reference at (1) is out of bounds (2 > 1) in dimension 1
./RADAU-INSTRUCOES.for:230:25:

       f(2)=-bb*x(2)+ff*x(3)-x(2)
                         1
Warning: Array reference at (1) is out of bounds (3 > 1) in dimension 1
./RADAU-INSTRUCOES.for:230:17:

       f(2)=-bb*x(2)+ff*x(3)-x(2)
                 1
Warning: Array reference at (1) is out of bounds (2 > 1) in dimension 1
./RADAU-INSTRUCOES.for:231:35:

       f(3)=5.d-4*bb*x(1)-ff*x(3)-x(2)
                                   1
Warning: Array reference at (1) is out of bounds (2 > 1) in dimension 1
./RADAU-INSTRUCOES.for:231:30:

       f(3)=5.d-4*bb*x(1)-ff*x(3)-x(2)
                              1
Warning: Array reference at (1) is out of bounds (3 > 1) in dimension 1
./RADAU-INSTRUCOES.for:285:29:

       WRITE(42,123)tm,X(1),X(2),X(3) ! each output_step (tm current) outputs the results
                             1
Warning: Array reference at (1) is out of bounds (2 > 1) in dimension 1
./RADAU-INSTRUCOES.for:285:34:

       WRITE(42,123)tm,X(1),X(2),X(3) ! each output_step (tm current) outputs the results
                                  1
Warning: Array reference at (1) is out of bounds (3 > 1) in dimension 1
./RADAU-INSTRUCOES.for:286:28:

       WRITE(*,123)tm,X(1),X(2),X(3)  ! save to file and display on screen
                            1
Warning: Array reference at (1) is out of bounds (2 > 1) in dimension 1
./RADAU-INSTRUCOES.for:286:33:

       WRITE(*,123)tm,X(1),X(2),X(3)  ! save to file and display on screen
                                 1
Warning: Array reference at (1) is out of bounds (3 > 1) in dimension 1
ra27.for:131:72:

      +OUTPUT(NF,NS,X,V,F1,T,TM,TF,.not.nsf,.FALSE.,iemax)
                                                                        1
Warning: More actual than formal arguments in procedure call at (1)
ra27.for:138:72:

       call OUTPUT(NF,NS,X,V,F1,T,TM,TF,.not.nsf,.FALSE.,iemax)
                                                                        1
Warning: More actual than formal arguments in procedure call at (1)
ra27.for:264:72:

       CALL OUTPUT(NF,NS,X,V,F1,T,TM,TF,.not.nsf,NPER,iemax)
                                                                        1
Warning: More actual than formal arguments in procedure call at (1)
/tmp/cc9kC96U.o: Na função "MAIN__":
RADAU-INSTRUCOES.for:(.text+0x36e9): referência não definida para "ra27_"
collect2: error: ld returned 1 exit status



reading the errors, you can see that one of the problems is the dimension of the variable x, which was defined in x (1).
However after defining the variable x as x (3), the following errors appear in ra27.for:

Code:
ra27.for:107:18:

       CALL FORCE (X, V, ZERO, F1)
                  1
Warning: Actual argument contains too few elements for dummy argument ‘x’ (1/3) at (1) [-Wargument-mismatch]
ra27.for:131:19:

      +OUTPUT(NF,NS,X,V,F1,T,TM,TF,.not.nsf,.FALSE.,iemax)
                   1
Warning: Actual argument contains too few elements for dummy argument ‘x’ (1/3) at (1) [-Wargument-mismatch]
ra27.for:138:24:

       call OUTPUT(NF,NS,X,V,F1,T,TM,TF,.not.nsf,.FALSE.,iemax)
                        1
Warning: Actual argument contains too few elements for dummy argument ‘x’ (1/3) at (1) [-Wargument-mismatch]
ra27.for:264:24:

       CALL OUTPUT(NF,NS,X,V,F1,T,TM,TF,.not.nsf,NPER,iemax)
                        1
Warning: Actual argument contains too few elements for dummy argument ‘x’ (1/3) at (1) [-Wargument-mismatch]
ra27.for:268:18:

   78  CALL FORCE (X,V,TM,F1)
                  1
Warning: Actual argument contains too few elements for dummy argument ‘x’ (1/3) at (1) [-Wargument-mismatch]
/tmp/ccfQadBO.o: Na função "MAIN__":
RADAU-INSTRUCOES.for:(text+0x36e9): referência não definida para "ra27_"
collect2: error: ld returned 1 exit status

RADAU-INSTRUCOES.FOR (WITHOUT CHANGES IN X):
Code:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

c: here is this trivial prog that uses the package: RA27.for (next to
c: Ra15 which is in fact an integration routine that uses here)
c: Still this black box is modified slightly
c: View original: Everhart article data output.
c: User must give Force.for and as an example goes the OUTPUT.for (routine
c: data output)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INCLUDE 'ra27.for'  ! include the package RA27.for

C: Just copy the RA27.for file to the same folder you are working on (same as this main program)
      
      IMPLICIT REAL*8(A-H,O-Z)

      logical fixed       
      CHARACTER FName1*30
      common/ctes/aa,bb,cc,dd,ee,ff
      data aa,bb,cc,dd,ee,ff/1.d0,0.2d0,0.1d0,0.6d0,-2.d0,-1.d0/

     
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                          C 
C                                                                          C
C                  C
C                                                 
C                                                                          C                                
C                             C                                            
C                                                                          C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C    dimension x: coordinate and v: velocity
c: in its case x = vector (x, y, z)
c: v = vector velocity of x, y, z: you do not need, your units are from 1st order

      DIMENSION X(3),v(3)   !  Assuming 3 first-order equations
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C NV 
C     
C      
C                                             
c: Nv: number of simultaneous events (maximum is 18, but
c: user can move (F1, FJ)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
     
      NV=3
      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C NCLASS      describes the type of system to be integrated 
C               y'=F(y,t)   NCLASS=1
C               y"=F(y,t)   NCLASS= -2
C               y"=F(y',y,t)  NCLASS=2
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      
      NCLASS=1  ! your system is the first order

c: eLagrange equations always sound first order, so nclass = 1

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C NOR  order of the integrator ie (15 or 27)
C  In case here, always use NOR = 15 (Radau15)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      
      NOR=15                                                              
      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C LL  controls to a certain extent accuracy (sequence size: see article). 
C       SS=10**(-LL) .
C      A good value of LL e ': half of NOR (maximum I think it is 12).
C    LL large-> high precision

C      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      
      LL=10 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C XL  if LL <= 0 then XL forms a constant-length sequence
C The smallest XL toleravel: XL = 0.01
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      XL=0.01D0

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C TF  integration interval TF = t (final) - t (initial)
C Logic that may be negative: integrate to Back             
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      
      tm=0.d0  ! time starts 
      TF=10.D0   ! final integration time

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     
C    X e V : position and initial velocity given by the user
C 
C
C     Entre aqui X(t=To) e V(t=To)=X' 
c     In case of first order, you do not need to give V (t = to)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C PI   with 32 digits. (German things)
C
      PI=3.1415926535897932384626433832795D0
C:INICIAL CONDITIONS of the variables x (1) x (2), x (3) (speed you do not have, do not need)

	x(1)=0.2D0  ! would be x (1) at t = t0
	x(2)=0.3d0  ! x(2) em t0
	x(3)=-0.05d0 ! x(3) em t0

c:	  ! in your case you do not need speeds

	    
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C OUTPUT_STEP: tells the integrator: exit result. to each (APROX) OUTPUT_step (time t)
C OUTPUT and called OUTPUT_STEP
C
C fixed: true if I want to exit in exactly OUTPUT_STEP
C
C: false if I do not ask too much, this is completed
c: sequence releases an output
C: Whenever possible use FIXED = .false. (much faster)

C		  
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                                                         
      OUTPUT_STEP=1.5D0  ! 


c: in fact it will not be exactly 1.5 (eh approx), if you want exactly 1.5, then change fixed to .true,
c: but alas, it gets a little slower, of course.
      
c:      fixed=.true.  ! here every exact output_step output the results, fixed pitch  
       fixed=.false.  ! in this case each output_step output the results (approximately)
c: use fixed = false, the integrator walks much faster

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC      
C
C  Filename for Output file: from the file where your outputs will be saved
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c:      FName1='/tmp/bla2.xyz'
         
        FName1='apg.dat'  !name of the file, to save the outputs                                                  
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      open (42,file=FName1, status='unknown') 

      write(42,*) ' Ordem da Integracao = ',NOR,
     +' SequenzSize = ', LL
      
      write(*,*) ' Ordem da Integracao = ',NOR,
     +' LL (precisao) = ', LL
    
      
 123  FORMAT(7F20.10)
   5  FORMAT (A)

      WRITE(*,5) ' Favor esperar: estou calculando ........'

      IF (NOR.eq.15)
     + CALL RA15 (X,V,TF,XL,LL,NV,NCLASS,NOR,OUTPUT_STEP,fixed)
      
      IF (NOR.gt.15)
     + CALL RA27 (X,V,TF,XL,LL,NV,NCLASS,NOR,OUTPUT_STEP,fixed)
      
    
      Close(42)

      write(*,5) ' Pronto, FIM '
      
      STOP
      END
c:ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C
C    C Now provide the Force.for subroutine, where it must conform to
C Nclass that you have defined before.
c
c
C SUBROUTINE FORCE (X, V, TM, F)
C X, V: see previous definition
C TM: Current time
C F: forca: vector of NV components
c: in general, Lagrange has no time = tm in the second limb
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE FORCE(X,V,TM,F)
      
          IMPLICIT REAL*8(A-H,O-Z)
          DIMENSION X(1),V(1),F(3)
	    common/ctes/aa,bb,cc,dd,ee,ff

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C     
C     Here you should give the components of the velocities of vector X (be careful
c: there must be agreement with NCLASS from the Main prog).
c:
c: F (1) = f1 (x, t)
c: F (2) = f2 (x, t)
c: ..............
c: F (nv) = fnv (x, t)

c Let X = [x (1), x (2), x (3)] be the easiest thing to do here:
c (1) = dx (1) / dt, f (2) = dx (2) / dt and f (3) = dx (3 / dt:
c: aa, bb, ff are constants given in prog. main

      f(1)=-aa*x(1)-bb*x(2)-0.002d0*x(3) ! f (i) are the first derivatives of x, y, z
      f(2)=-bb*x(2)+ff*x(3)-x(2)
      f(3)=5.d-4*bb*x(1)-ff*x(3)-x(2)
   
      RETURN
      END        
 
C ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc    
C     
C This is the routine that prints the outputs on the screen and in the file apg.apg
c
C			 SUBROUTINE OUTPUT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
      SUBROUTINE OUTPUT(NF,NS,X,V,F,T,TM,TF,FIRST,LAST)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C NF   :Number of iterations made by Radau 
C
C NS   : number of steps in time 
C
C X,V   X and V : said before
C
C  F   power vector
C
C  T   size of the next step (time t)
C
C TM   current time
C
C TF   TF=t(final) - t (initial) may be negative
C     
C FIRST logical variable: truth only at the beginning
C The integrator can test several times at the beginning until it finds a good
c: to start
c:  
C      
cc All this information can be ignored in a first read, Just remember tm = current time         
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                                                     
      IMPLICIT REAL*8(A-H,O-Z)
      INTEGER NF,NS,N
      LOGICAL LAST,FIRST
 
      DIMENSION X(1),V(1),F(1)  ! as this is a fortran 77, in general the vectors of
c: subroutines, are scaled to 1 only. If you are going to use Gfortran, then you should set the maximum, like
c: was defined in the main programl (3)
         

   5  FORMAT (A)
 123  FORMAT(F16.8, 3d21.10)  
      
C  Here I am writing the outputs: screen and file in unit 42 (attention
c: this 42 must agree with the OPEN: prog. Main)

C: if there is an integral integral.teste a preciso com several LL !!!             
      WRITE(42,123)tm,X(1),X(2),X(3) ! each output_step (tm current) outputs the results
      WRITE(*,123)tm,X(1),X(2),X(3)  ! save to file and display on screen
  
      RETURN
  
      END        
      
cccccccccccccccccccccccccccccccccccccccccccccccccccc Fim cccccccccccccccccccccccccccccccccccc


ra27.for:
Code:
C  $$$$  RA15T.FOR
C
      SUBROUTINE RA15(X,V,TF,XL,LL,NV,NCLASS,NOR,OS,FIXED)
C  Integrator RADAU by E. Everhart, Physics Department, University of Denver
C  This 15th-order version, called RA15, is written out for faster execution.
C  y'=F(y,t) is  NCLASS=1,  y"=F(y,t) is NCLASS= -2,  y"=F(y',y,t) is NCLASS=2
C  TF is t(final) - t(initial). It may be negative for backward integration.
C  NV = the number of simultaneous differential equations.
C  The dimensioning below assumes NV will not be larger than 18.
C  LL controls sequence size. Thus SS=10**(-LL) controls the size of a term.
C  A typical LL-value is in the range 6 to 12 for this order 11 program.
C  However, if LL.LT.0 then XL is the constant sequence size used.
C  X and V enter as the starting position-velocity vector, and are output as
C  the final position-velocity vector.
C  Integration is in double precision. A 64-bit double-word is assumed.
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL *4 TVAL,PW
      REAL*8 NEXTOUTPUT
      DIMENSION X(1),V(1),F1(18),FJ(18),C(21),D(21),R(21),Y(18),Z(18),
     A     B(7,18),G(7,18),E(7,18),BD(7,18),H(8),W(7),U(7),NW(8)
      
c: 25/11/98 introduzi iemax: introduce iemax: stop integration if eccentricity reaches 0.7
c: iemax goes into the OUTPUT_STEP subroutine
      
      LOGICAL NPQ,NSF,NPER,NCL,NES,fixed
      DATA NW/0,0,1,3,6,10,15,21/
      DATA ZERO, HALF, ONE,SR/0.0D0, 0.5D0, 1.0D0,1.4D0/
C  These H values are the Gauss-Radau spacings, scaled to the range 0 to 1,
C  for integrating to order 15.
      DATA H/         0.D0, .05626256053692215D0, .18024069173689236D0,
     A.35262471711316964D0, .54715362633055538D0, .73421017721541053D0,
     B.88532094683909577D0, .97752061356128750D0/
C  The sum of the H-values should be 3.73333333333333333

c: 26/11/98 : large eccentricity control with IEMAX
      iemax=0
       
      NPER=.FALSE.
      NSF=.FALSE.
      NCL=NCLASS.EQ.1
      NPQ=NCLASS.LT.2
      Out=OS
C  y'=F(y,t)  NCL=.TRUE.    y"=F(y,t)  NCL=.FALSE.   y"=F(y',y,t) NCL=.FALSE.
C  NCLASS=1   NPQ=.TRUE.    NCLASS= -2 NPQ=.TRUE.    NCLASS= 2    NPQ=.FALSE.
C  NSF is .FALSE. on starting sequence, otherwise .TRUE.
C  NPER is .TRUE. only on last sequence of the integration.
C  NES is .TRUE. only if LL is negative. Then the sequence size is XL.
      DIR=ONE
      IF(TF.LT.ZERO) DIR=-ONE
      NES=LL.LT.0
      XL=DIR*DABS(XL)
      PW=1./9.
C  Evaluate the constants in the W-, U-, C-, D-, and R-vectors
      DO 14 N=2,8
      WW=N+N*N
      IF(NCL) WW=N
      W(N-1)=ONE/WW
      WW=N
  14  U(N-1)=ONE/WW
      DO 22 K=1,NV
      IF(NCL) V(K)=ZERO
      DO 22 L=1,7
      BD(L,K)=ZERO
  22  B(L,K)=ZERO
      W1=HALF
      IF(NCL) W1=ONE
      C(1)=-H(2)
      D(1)=H(2)
      R(1)=ONE/(H(3)-H(2))
      LA=1
      LC=1
      DO 73 K=3,7
      LB=LA
      LA=LC+1
      LC=NW(K+1)
      C(LA)=-H(K)*C(LB)
      C(LC)=C(LA-1)-H(K)
      D(LA)=H(2)*D(LB)
      D(LC)=-C(LC)
      R(LA)=ONE/(H(K+1)-H(2))
      R(LC)=ONE/(H(K+1)-H(K))
      IF(K.EQ.3) GO TO 73
      DO 72 L=4,K
      LD=LA+L-3
      LE=LB+L-4
      C(LD)=C(LE)-H(K)*C(LE+1)
      D(LD)=D(LE)+H(L-1)*D(LE+1)
  72  R(LD)=ONE/(H(K+1)-H(L-1))
  73  CONTINUE
      SS=10.**(-LL)
C  The statements above are used only once in an integration to set up the
C  constants. They use less than a second of execution time.  Next set in
C  a reasonable estimate to TP based on experience. Same sign as DIR.
C  An initial first sequence size can be set with XL even with LL positive.
      TP=0.1D0*DIR
      IF(XL.NE.ZERO) TP=XL
      IF(NES) TP=XL
      IF(TP/TF.GT.HALF) TP=HALF*TF
      NCOUNT=0

C  An * is the symbol for writing on the monitor. The printer is unit 4.
C  Line 4000 is the starting place of the first sequence.
4000  NS=0
      NF=0
      NI=6
      TM=ZERO
      CALL FORCE (X, V, ZERO, F1)
      NF=NF+1
C Line 722 is begins every sequence after the first. First find new beta-
C  values from the predicted B-values, following Eq. (2.7) in text.
 722  DO 58 K=1,NV
      G(1,K)=B(1,K)+D(1)*B(2,K)+
     X  D(2)*B(3,K)+D(4)*B(4,K)+D( 7)*B(5,K)+D(11)*B(6,K)+D(16)*B(7,K)
      G(2,K)=            B(2,K)+
     X  D(3)*B(3,K)+D(5)*B(4,K)+D( 8)*B(5,K)+D(12)*B(6,K)+D(17)*B(7,K)
      G(3,K)=B(3,K)+D(6)*B(4,K)+D( 9)*B(5,K)+D(13)*B(6,K)+D(18)*B(7,K)
      G(4,K)=            B(4,K)+D(10)*B(5,K)+D(14)*B(6,K)+D(19)*B(7,K)
      G(5,K)=                         B(5,K)+D(15)*B(6,K)+D(20)*B(7,K)
      G(6,K)=                                      B(6,K)+D(21)*B(7,K)
  58  G(7,K)=                                                   B(7,K)
      T=TP
      T2=T*T
      IF(NCL) T2=T
      TVAL=DABS(T)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c: 25/11/98  
c: here are the OUTPUT_STEP calls: Then, as in this BERND3, I interrupt the integration
c: if eccentricity reaches 0.7 (in this case iemax = 222), make RA15 return 

      IF(fixed.and.DABS(DIR*TM-OUT+OS).le.1.d-8) call  
     +OUTPUT(NF,NS,X,V,F1,T,TM,TF,.not.nsf,.FALSE.,iemax)
           
      if(iemax.eq.222) return
c: if iemax = 222, I can stop current integration and restart new orbit.
           
      IF (fixed.or.(out-os-dir*tm).gt.1.d-8)
     +go to 246
      call OUTPUT(NF,NS,X,V,F1,T,TM,TF,.not.nsf,.FALSE.,iemax)
      Out=out+os 
      
      if(iemax.eq.222) return
c: if iemax = 222, stop the current integration and start new orbit
      
  246 continue
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
      
C  Loop 175 is 6 iterations on first sequence and two iterations therafter.
      DO 175 M=1,NI
C  Loop 174 is for each substep within a sequence.
      DO 174 J=2,8
      JD=J-1
      JDM=J-2
      S=H(J)
      Q=S
      IF(NCL) Q=ONE
C  Use Eqs. (2.9) and (2.10) of text to predict positions at each aubstep.
C  These collapsed series are broken into two parts because an otherwise
C  excellent  compiler could not handle the complicated expression.
      DO 130 K=1,NV
      A=W(3)*B(3,K)+S*(W(4)*B(4,K)+S*(W(5)*B(5,K)+S*(W(6)*B(6,K)+
     V   S*W(7)*B(7,K))))
      Y(K)=X(K)+Q*(T*V(K)+T2*S*(F1(K)*W1+S*(W(1)*B(1,K)+S*(W(2)*B(2,K)
     X  +S*A))))
      IF(NPQ) GO TO 130
C  Next are calculated the velocity predictors need for general class II.
      A=U(3)*B(3,K)+S*(U(4)*B(4,K)+S*(U(5)*B(5,K)+S*(U(6)*B(6,K)+
     T    S*U(7)*B(7,K))))
      Z(K)=V(K)+S*T*(F1(K)+S*(U(1)*B(1,K)+S*(U(2)*B(2,K)+S*A)))
 130  CONTINUE
C  Find forces at each substep.
      CALL FORCE(Y,Z,TM+S*T,FJ)
      NF=NF+1
      DO 171 K=1,NV
C  Find G-value for the force FJ found at the current substep. This
C  section, including the many-branched GOTO, uses Eq. (2.4) of text.
      TEMP=G(JD,K)
      GK=(FJ(K)-F1(K))/S
      GO TO (102,102,103,104,105,106,107,108),J
 102  G(1,K)=GK
      GO TO 160
 103  G(2,K)=(GK-G(1,K))*R(1)
      GO TO 160
 104  G(3,K)=((GK-G(1,K))*R(2)-G(2,K))*R(3)
      GO TO 160
 105  G(4,K)=(((GK-G(1,K))*R(4)-G(2,K))*R(5)-G(3,K))*R(6)
      GO TO 160
 106  G(5,K)=((((GK-G(1,K))*R(7)-G(2,K))*R(8)-G(3,K))*R(9)-
     X     G(4,K))*R(10)
      GO TO 160
 107  G(6,K)=(((((GK-G(1,K))*R(11)-G(2,K))*R(12)-G(3,K))*R(13)-
     X     G(4,K))*R(14)-G(5,K))*R(15)
      GO TO 160
 108  G(7,K)=((((((GK-G(1,K))*R(16)-G(2,K))*R(17)-G(3,K))*R(18)-
     X     G(4,K))*R(19)-G(5,K))*R(20)-G(6,K))*R(21)
C  Upgrade all B-values.
 160  TEMP=G(JD,K)-TEMP
      B(JD,K)=B(JD,K)+TEMP
C  TEMP is now the improvement on G(JD,K) over its former value.
C  Now we upgrade the B-value using this dfference in the one term.
C  This section is based on Eq. (2.5).
      GO TO (171,171,203,204,205,206,207,208),J
 203  B(1,K)=B(1,K)+C(1)*TEMP
      GO TO 171
 204  B(1,K)=B(1,K)+C(2)*TEMP
      B(2,K)=B(2,K)+C(3)*TEMP
      GO TO 171
 205  B(1,K)=B(1,K)+C(4)*TEMP
      B(2,K)=B(2,K)+C(5)*TEMP
      B(3,K)=B(3,K)+C(6)*TEMP
      GO TO 171
 206  B(1,K)=B(1,K)+C(7)*TEMP
      B(2,K)=B(2,K)+C(8)*TEMP
      B(3,K)=B(3,K)+C(9)*TEMP
      B(4,K)=B(4,K)+C(10)*TEMP
      GO TO 171
 207  B(1,K)=B(1,K)+C(11)*TEMP
      B(2,K)=B(2,K)+C(12)*TEMP
      B(3,K)=B(3,K)+C(13)*TEMP
      B(4,K)=B(4,K)+C(14)*TEMP
      B(5,K)=B(5,K)+C(15)*TEMP
      GO TO 171
 208  B(1,K)=B(1,K)+C(16)*TEMP
      B(2,K)=B(2,K)+C(17)*TEMP
      B(3,K)=B(3,K)+C(18)*TEMP
      B(4,K)=B(4,K)+C(19)*TEMP
      B(5,K)=B(5,K)+C(20)*TEMP
      B(6,K)=B(6,K)+C(21)*TEMP
 171  CONTINUE
 174  CONTINUE
      IF(NES.OR.M.LT.NI) GO TO 175
C  Integration of sequence is over. Next is sequence size control.
      HV=ZERO
      DO 635 K=1,NV
 635  HV=DMAX1(HV,DABS(B(7,K)))
      HV=HV*W(7)/TVAL**7
 175  CONTINUE
      IF (NSF) GO TO 180
      IF(.NOT.NES) TP=(SS/HV)**PW*DIR
      IF(NES) TP=XL
      IF(NES) GO TO 170
      IF(TP/T.GT.ONE) GO TO 170
   8  FORMAT (A,2X,2I2,2D18.10)
      TP=.8D0*TP
      NCOUNT=NCOUNT+1
      IF(NCOUNT.GT.10) RETURN
     
      
C  Restart with 0.8x sequence size if new size called for is smaller than
C  originally chosen starting sequence size on first sequence.
      GO TO 4000
 170  NSF=.TRUE.
C Loop 35 finds new X and V values at end of sequence using Eqs. (2.11),(2.12)
 180  DO 35 K=1,NV
      X(K)=X(K)+V(K)*T+T2*(F1(K)*W1+B(1,K)*W(1)+B(2,K)*W(2)+B(3,K)*W(3)
     X    +B(4,K)*W(4)+B(5,K)*W(5)+B(6,K)*W(6)+B(7,K)*W(7))
      IF(NCL) GO TO 35
      V(K)=V(K)+T*(F1(K)+B(1,K)*U(1)+B(2,K)*U(2)+B(3,K)*U(3)
     V    +B(4,K)*U(4)+B(5,K)*U(5)+B(6,K)*U(6)+B(7,K)*U(7))
  35  CONTINUE
      TM=TM+T
      NS=NS+1
C  Return if done.
      IF(.NOT.NPER) GO TO 78
      CALL OUTPUT(NF,NS,X,V,F1,T,TM,TF,.not.nsf,NPER,iemax)
      RETURN
C  Control on size of next sequence and adjust last sequence to exactly
C  cover the integration span. NPER=.TRUE. set on last sequence.
  78  CALL FORCE (X,V,TM,F1)
      NF=NF+1
      IF(NES) GO TO 341
      TP=DIR*(SS/HV)**PW
      IF(TP/T.GT.SR) TP=T*SR
 341  IF(NES) TP=XL
      IF(DIR*(TM+TP).LT.DIR*TF-1.D-8) GO TO 66
      TP=TF-TM
      NPER=.TRUE.

 66   IF (.not.fixed.or.DIR*(TM+TP).LT.OUT-1.D-8) GO TO 77
      TP=DIR*OUT-TM
      OUT=OUT+OS 
      
C  Now predict B-values for next step. The predicted values from the preceding
C  sequence were saved in the E-matrix. Te correction BD between the actual
C  B-values found and these predicted values is applied in advance to the
C  next sequence. The gain in accuracy is significant. Using Eqs. (2.13):
  77  Q=TP/T
      DO 39 K=1,NV
      IF(NS.EQ.1) GO TO 31
      DO 20 J=1,7
  20  BD(J,K)=B(J,K)-E(J,K)
  31  E(1,K)=      Q*(B(1,K)+ 2.D0*B(2,K)+ 3.D0*B(3,K)+
     X           4.D0*B(4,K)+ 5.D0*B(5,K)+ 6.D0*B(6,K)+ 7.D0*B(7,K))
      E(2,K)=                Q**2*(B(2,K)+ 3.D0*B(3,K)+
     Y           6.D0*B(4,K)+10.D0*B(5,K)+15.D0*B(6,K)+21.D0*B(7,K))
      E(3,K)=                             Q**3*(B(3,K)+
     Z           4.D0*B(4,K)+10.D0*B(5,K)+20.D0*B(6,K)+35.D0*B(7,K))
      E(4,K)=   Q**4*(B(4,K)+ 5.D0*B(5,K)+15.D0*B(6,K)+35.D0*B(7,K))
      E(5,K)=                Q**5*(B(5,K)+ 6.D0*B(6,K)+21.D0*B(7,K))
      E(6,K)=                             Q**6*(B(6,K)+ 7.D0*B(7,K))
      E(7,K)=                                           Q**7*B(7,K)
      DO 39 L=1,7
  39  B(L,K)=E(L,K)+BD(L,K)
C  Two iterations for every sequence after the first.
      NI=2
      GO TO 722
      END
.
 
Try the following:
Those arrays that were declared with full length, x(3), in the main program but with just on, x(1), in the subroutines, declared them with full length in every subroutine, too.
 
You could also define the arrays in a common block used by all modules. Then it will use the same memory and format for all your arrays

Bill
Lead Application Developer
New York State, USA
 
Status
Not open for further replies.

Part and Inventory Search

Sponsor

Back
Top