COMPILAR PROGRAMA FORTRAN RADAU

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 

albfneto
(usa openSUSE)

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.

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),"