COMPILAR PROGRAMA FORTRAN RADAU

1. COMPILAR PROGRAMA FORTRAN RADAU

Victor Hugo Mota
pistoler

(usa Ubuntu)

Enviado em 03/10/2018 - 17:45h

Oi pessoal, eu não to conseguindo compilar um programa aqui (RADAU-INSTRUCOES.for), que usa o programa radau ra15, (ra27.for)
Quando compilo o programa RADAU-INSTRUCOES.for, aparece os seguintes erros;
./RADAU-INSTRUCOES.for:228:38:

f(1)=-aa*x(1)-bb*x(2)-0.002d0*x(3) ! os f(i) sao as derivadas primeiras de x,y,z
1
Warning: Array reference at (1) is out of bounds (3 > 1) in dimension 1
./RADAU-INSTRUCOES.for:228:25:

f(1)=-aa*x(1)-bb*x(2)-0.002d0*x(3) ! os f(i) sao as derivadas primeiras de x,y,z
1
Warning: Array reference at (1) is out of bounds (2 > 1) in dimension 1
./RADAU-INSTRUCOES.for:229: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:229: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:229: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:230: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:230: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:287:29:

WRITE(42,123)tm,X(1),X(2),X(3) ! a cada output_step (tm atual) saem os results
1
Warning: Array reference at (1) is out of bounds (2 > 1) in dimension 1
./RADAU-INSTRUCOES.for:287:34:

WRITE(42,123)tm,X(1),X(2),X(3) ! a cada output_step (tm atual) saem os results
1
Warning: Array reference at (1) is out of bounds (3 > 1) in dimension 1
./RADAU-INSTRUCOES.for:288:28:

WRITE(*,123)tm,X(1),X(2),X(3) ! gravo no arquivo e exibo na tela
1
Warning: Array reference at (1) is out of bounds (2 > 1) in dimension 1
./RADAU-INSTRUCOES.for:288:33:

WRITE(*,123)tm,X(1),X(2),X(3) ! gravo no arquivo e exibo na tela
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/ccinIWAr.o: Na função "MAIN__":
RADAU-INSTRUCOES.for:(.text+0x36e9): referência não definida para "ra27_"
collect2: error: ld returned 1 exit status

lendo os erros , dá para perceber que um dos problemas é a dimensão da variavel x, que estava definida em x(1). Porém depois de definir a variavel x como x(3), aparece os seguintes erros no ra27.for:
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/cczaNxX4.o: Na função "MAIN__":
RADAU-INSTRUCOES.for:(.text+0x36e9): referência não definida para "ra27_"
collect2: error: ld returned 1 exit status

Alguém pode me explicar o que está de errado?


RADAU-INSTRUCOES.FOR (SEM ALTERAÇÕES EM X)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

c: aqui esta um prog. trivial que usa o pacotao: RA27.for (contem junto o
c: Ra15 que eh de fato a rotina de integracao que uso aqui)
c: Mesmo assim esta caixa preta esta modificada um pouco
c: Ver original : artigo do Everhart saida de dados.
c: Usuario deve dar Force.for e como exemplo vai o OUTPUT.for (rotina
c: de saida de dados)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

INCLUDE 'ra27.for' ! inclua o pacote RA27.for

C: Apenas copie o arquivo RA27.for na mesma pasta que voce esta trabalhando( a mesma deste prog. principal)

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 Aqui : algumas coisas que o usuario deve dar C
C
C C
C C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C dimensione o x: coordenada e v: velocidade
c: no seu caso x= vetor (x,y,z)
c: v= vetor velocidade de x,y,z : voce nao precisa, suas eqcoes sao de 1-meira ordem

DIMENSION X(3),v(3) ! Supondo 3 eqcoes de primeira ordem
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C NV
C
C
C
c: Nv: numero de eqcoes simultanaes (maximo esta em 18, mas
c: usuario pode mexer (F1, FJ )
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

NV=3

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C NCLASS descreve o tipo de sistema a ser integrado
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 ! seu sistema eh de primeira ordem

c: eqcoes de Lagrange soa sempre de primeira ordem, entao nclass=1

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C NOR ordem do integrador ou seja (15 ou 27)
C No caso aqui, use sempre NOR=15 (Radau15)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

NOR=15

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C LL controla de certo modo a precisao (tamanho da sequencia: ver artigo).
C SS=10**(-LL) .
C Um bom valor de LL e': metade de NOR (maximo acho que é 12).
C LL grande-> alta precisao

C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

LL=10


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C XL se LL<=0 entao XL forma uma sequencia de tamanho constante
C O menor XL toleravel: XL = 0.01
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

XL=0.01D0

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C TF intervalo de integracao TF=t(final) - t(inicial)
C Logico que pode ser negativo: integrar p/Tras
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

tm=0.d0 ! tempo inicial
TF=10.D0 ! tempo final de integracao

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
C X e V : posicao e velocidade inicial dadas pelo usuario
C
C
C Entre aqui X(t=To) e V(t=To)=X'
c No caso de primeira ordem, voce nao precisa dar V(t=to)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C PI com 32 digitos. (coisas de alemao)
C
PI=3.1415926535897932384626433832795D0
C: CONDICOES INICIAIS das variaveis x(1) x(2),x(3) (velocidade voce nao tem, nao precisa)

x(1)=0.2D0 ! seria o x(1) em t=t0
x(2)=0.3d0 ! x(2) em t0
x(3)=-0.05d0 ! x(3) em t0

c: ! no seu caso voce nao precisa de velocidades

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C OUTPUT_STEP:diz ao integrador: saia result. a cada (APROX) OUTPUT_step (tempo t)
C OUTPUT e chamado acada OUTPUT_STEP
C
C fixed : verdade se quero saida em exatamente OUTPUT_STEP
C
C : falso se nao faco muita questao, isto e'qdo completada uma
c: sequencia libero uma saida
C : Sempre que possivel use FIXED=.false. (muito mais rapido)
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

OUTPUT_STEP=1.5D0 ! significa: a cada 1.5 (anos, dias, ou horas) saem os x(i)


c: na verdade nao sera exatamente 1.5 (eh aproximado), se se quiser exatamente 1.5, entao mude fixed para .true.,
c: mas ai, fica um pouco mais lento, obvio.

c: fixed=.true. ! aqui a cada exatos output_step saem os resultados , passo fixo
fixed=.false. ! neste caso a cada output_step saem os resultados (aproximadamente)
c: usar fixed=false, o integrador caminha muito mais rapidamente

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Filename para Output file: de o arquivo onde serao guardados sua saidas
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c: FName1='/tmp/bla2.xyz'

FName1='apg.dat' ! nome do arquivo , p/ guardar as saidas
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 Agora forneca subrotina Force.for, onde ela deve estar de acordo com
C Nclass que voce definiu antes.
C
C
C SUBROUTINE FORCE(X,V,TM,F)
C X,V : ver definicao anterior
C TM : tempo atual
C F : forca : vetor de NV componentes
c: em geral eqcoes de Lagrange nao tem tempo= tm no segundo membro

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 Aqui voce deve dar as componentes das velocidades do vetor X (cuidado que
c: deve haver acordo com NCLASS do prog Principal).
c:
c: F(1)= f1(x,t)
c: F(2)= f2(x,t)
c: ..............
c: F(nv)=fnv(x,t)

c: P/ facilitar, vou fazer aqui o dicionario das variaveis : Seja X=[x(1),x(2),x(3) ] entao
c: f(1)= dx(1)/dt, f(2)=dx(2)/dt e f(3)= dx(3/dt : eqcoes de Lagrange ficticias
c: aa,bb,ff sao constantes dadas no prog. principal

f(1)=-aa*x(1)-bb*x(2)-0.002d0*x(3) ! os f(i) sao as derivadas primeiras de 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 Esta eh a rotina que imprime oas saidas na tela e no arquivo apg.apg
C
C SUBROUTINE OUTPUT
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
SUBROUTINE OUTPUT(NF,NS,X,V,F,T,TM,TF,FIRST,LAST)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C NF :numero iteracoes feito pelo Radau
C
C NS : numero de passos no tempo
C
C X,V X and V : dito antes
C
C F vetor forca
C
C T tamanho do proximo passo (tempo t)
C
C TM tempo atual
C
C TF TF=t(final) - t(inicial) pode ser negativo
C
C FIRST variavel logica: verdade so no inicio
C O integrador pode testar varias vezes no inicio ate achar um bom
c: p/ comecar
C
C LAST : logica : verdade so no fim do ultimo passo.
c:
C
C Todas estas informacoes, podem ser ignoradas numa primeira lida, Apenas lembrar tm=tempo atual
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

IMPLICIT REAL*8(A-H,O-Z)
INTEGER NF,NS,N
LOGICAL LAST,FIRST

DIMENSION X(1),V(1),F(1) ! como isto e'um fortran 77, em geral os vetores de
c: sub rotinas, sao dimensionados em 1 apenas. Se voce for usar Gfortran, entao deve botar dimensao maxima, como
c: foi definida no prog princiapl (3)


5 FORMAT (A)
123 FORMAT(F16.8, 3d21.10)

C Aqui estou escrevendo as saidas: tela e arquivo na unidade 42 (atencao
c: este 42 deve concordar com o OPEN: prog. Principal)

C: se houver uma integral primeira.teste a precisao com varios LL!!!
WRITE(42,123)tm,X(1),X(2),X(3) ! a cada output_step (tm atual) saem os results
WRITE(*,123)tm,X(1),X(2),X(3) ! gravo no arquivo e exibo na tela

RETURN

END

cccccccccccccccccccccccccccccccccccccccccccccccccccc Fim cccccccccccccccccccccccccccccccccccc




ra27.for;
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: parar integracao se excentricidade chegar a 0.7
c: iemax vai na sub-rotina OUTPUT_STEP

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 : controle de excentricidade grande com 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: aqui existem as chamadas do OUTPUT_STEP: Entao como neste BERND3 interrompo a integracao
c: se excentricidade chega a 0.7 (neste caso iemax=222), faco RA15 retornar .

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: se iemax=222 ,posso parar integracao atual e recomecar nova orbita

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: se iemax=222, parar a atual integracao e comecar nova orbita

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



  


2. Re: COMPILAR PROGRAMA FORTRAN RADAU

Alberto Federman Neto.
albfneto

(usa Sabayon)

Enviado em 03/10/2018 - 18:19h

Aquarde mais respostas de galera que conhece Fortran.
Não sou programador, nem profissional de TI,
¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨
Albfneto,
Ribeirão Preto, S.P., Brasil.
Usuário Linux, Linux Counter: #479903.
Distros Favoritas: Sabayon, Gentoo, openSUSE, Mageia e OpenMandriva.


3. Re: COMPILAR PROGRAMA FORTRAN RADAU

Rubens Antonio Condeles Júnior
Rcondeles

(usa openSUSE)

Enviado em 04/11/2018 - 10:36h

Olá,
O erro ocorre quando o programa principal chama a subrotina ra27. Em ra27.for esta subrotina não existe.
Comente as linhas com "!" (ou você pode excluí-las) :

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

agora o compilador irá gerar o programa executável.

Entretanto, ainda ocorrerão várias advertências ("Warning:").
Elas ocorrem pelo fato de que os vetores X, V e F foram declarados
com "shape" diferente nas subrotinas. No seu exemplo os vetores
X,V e F apresentam "shape" = 3

Para a correção faça as seguintes alterações:

Na subroutine force:

onde aparece "DIMENSION X(1),V(1),F(3)" mude para "DIMENSION X(3),V(3),F(3)"

Na subroutine output:

onde aparece "SUBROUTINE OUTPUT(NF,NS,X,V,F,T,TM,TF,FIRST,LAST)"
mude para "SUBROUTINE OUTPUT(NF,NS,X,V,F,T,TM,TF,FIRST,LAST,iemax)"

onde aparece "DIMENSION X(1),V(1),F(1)"
mude para "DIMENSION X(3),V(3),F(3)"


Na subroutine ra15:

onde aparece "DIMENSION X(1),V(1),F1(18),FJ(18),C(21),D(21),R(21),Y(18),Z(18),"
mude para "DIMENSION X(3),V(3),F1(18),FJ(18),C(21),D(21),R(21),Y(18),Z(18),"






Patrocínio

Site hospedado pelo provedor RedeHost.
Linux banner
Linux banner
Linux banner

Destaques

Artigos

Dicas

Tópicos

Top 10 do mês

Scripts