×
INTELLIGENT WORK FORUMS
FOR COMPUTER PROFESSIONALS

Are you a
Computer / IT professional?
Join Tek-Tips Forums!
• Talk With Other Members
• Be Notified Of Responses
• Keyword Search
Favorite Forums
• Automated Signatures
• Best Of All, It's Free!

*Tek-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

Jobs

 Forum Search FAQs Links MVPs

(OP)
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 --> fortran

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

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

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

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

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

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

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

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

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

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

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 --> fortran

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 --> fortran

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 --> fortran

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
.

RE: Compiler Fortran RADAU

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.

RE: Compiler Fortran RADAU

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
New York State, USA

Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

Red Flag Submitted

Thank you for helping keep Tek-Tips Forums free from inappropriate posts.
The Tek-Tips staff will check this out and take appropriate action.

Posting in the Tek-Tips forums is a member-only feature.

Click Here to join Tek-Tips and talk with other members! Already a Member? Login

Close Box

Join Tek-Tips® Today!

Join your peers on the Internet's largest technical computer professional community.
It's easy to join and it's free.

Here's Why Members Love Tek-Tips Forums:

• Talk To Other Members
• Notification Of Responses To Questions
• Favorite Forums One Click Access
• Keyword Search Of All Posts, And More...

Register now while it's still free!