NASA Ocean Color

ocssw V2022
       SUBROUTINE rk78(DER,T,TOUT,N,X,H,RELERR,ABSERR,SP)  
 C  VERSION OF 4/1/85
 C  PURPOSE
 C    RUNGE-KUTTA 7(8) METHOD AS GIVEN BY THE FEHLBERG COEFFICIENTS
 C    NASA TR R-287
 C  INPUT
 C    DER    = NAME OF SUBROUTINE CONTAINING THE DERIVATIVES
 C    T      = INITIAL INTEGRATION TIME
 C    TOUT   = FINAL INTEGRATION TIME
 C    N      = NUMBER OF EQUATIONS TO BE INTEGRATED
 C    X      = THE STATE TO BE INTEGRATED, DIMENSION 6 HERE 
 C    H      = INITIAL GUESS OF STEP-SIZE
 C    RELERR = RELATIVE TOLERANCE
 C    ABSERR = ABSOLUTE TOLERANCE
 C    SP     = APPROXIMATE TIME THAT USEROP WILL BE CALLED
 C  OUTPUT
 C    T      = AT TOUT, T=TOUT
 C    X      = STATE AT TOUT
 C  CALL SUBROUTINES
 C    DER, USEROP
 C  REFERENCES
 C    JPL EM 312/87-153, 20 APRIL 1987
 C    NASA TR R-287, OCTOBER 1968
 C  ANALYSIS
 C    J. H. KWOK - JPL
 C  PROGRAMMER
 C    J. H. KWOK - JPL
 C  PROGRAM MODIFICATIONS
 C
 C  DECLARED DER AS EXTERNAL TO ELIMINATE COMPILER WARNING.
 C  F. S. PATT, JANUARY 13, 1995.
 C
 C    NONE
 C  COMMENTS
 C    USER MUST CALL RK78CN BEFORE CALLING RK78.
 C    THE DIMENSIONS OF XDUM AND F CAN BE CHANGED FROM 6 TO ACCOMMODATE
 C    ANY SYSTEM OF DIFFERENTIAL EQUATIONS.  USEROP IS A USER OPTION
 C    ROUTINE WHERE THE HEALTH OF THE INTEGRATION CAN BE MONITORED BY
 C    SETTING SP=ZERO.  IT CAN ALSO BE USED TO CHECK CERTAIN CONDITIONS
 C    OF THE STATE, AND POSSIBLY RESTART WITH A NEW STATE OR STEP-SIZE
 C    BY SETTING ISTART=1 IN USEROP.
 C
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
       common/felcon/ch(13),alph(13),beta(13,12),ordrcp,norder,ntimes
       dimension xdum(6),f(6,13),x(*)
       EXTERNAL der
       DATA zero,one/0.d0,1.d0/
 C
 C *** FACT IS A SCALING FACTOR TO REDUCE THE ESTIMATED STEPSIZE TO AVOID
 C *** EXCEESIVE NUMBER OF REJECTION.
 C *** SMALL IS A SMALL NUMBER COMPARED TO THE SIZE OF T
 C
       DATA fact,small/0.7d0,1.d-8/
 C
 C *** INITIALIZATION
 C
       nstp=0
       nrej=0
       sdt=dsign(one,tout-t)   
       dt=dabs(h)*sdt
       spp=dabs(sp)*sdt
       tpr=t
       istart=0
       iflag=0
 C
 C *** CALL USEROP AT INITIAL T
 C
 C      CALL USEROP(T,X,N,DT,NREJ,NSTP,ISTART,IFLAG)
 C
 C *** START INTEGRATION AND SAVE INITIAL INFO
 C
   100 CONTINUE
       istart=0
       nstp=nstp+1
       tdum=t
       tpr=tpr+spp
       DO 110 i=1,n   
   110 xdum(i)=x(i)
 C
 C *** CHECK FOR T+DT PASSING TOUT
 C
       IF (dabs(dt).GT.dabs(tout-t)) dt=tout-t  
 C
 C *** CHECK FOR REACHING TOUT
 C
       IF (dabs(t-tout).LT.small) GO TO 900
 C
 C *** START FUNCTION EVALUATIONS
 C
       CALL der(t,x,f(1,1))
       DO 300 k=2,ntimes  
       kk=k-1
       DO 310 i=1,n   
       temp=zero
       DO 320 j=1,kk  
   320 temp=temp+beta(k,j)*f(i,j)
   310 x(i)=xdum(i)+dt*temp  
       t=tdum+alph(k)*dt 
       CALL der(t,x,f(1,k))  
   300 CONTINUE
 C
 C *** COMPUTE NEW STATE
 C
       DO 400 i=1,n  
       temp=zero 
       DO 410 l=1,ntimes 
   410 temp=temp+ch(l)*f(i,l)
   400 x(i)=xdum(i)+dt*temp
 C
 C *** START LOCAL TRUNCATION ERROR COMPUTATION
 C
       err=relerr
       DO 500 i=1,n  
       ter=dabs((f(i,1)+f(i,11)-f(i,12)-f(i,13))*ch(12)*dt)
       tol=dabs(x(i))*relerr+abserr
       const=ter/tol
       IF (const.GT.err) err=const
   500 CONTINUE  
 C
 C *** ESTIMATE NEW STEP-SIZE
 C
       dt=fact*dt*(one/err)**ordrcp
 C
 C *** IF IFLAG.EQ.1, ACCEPT THIS STEP UNCONDITIONALLY
 C
       IF (iflag.EQ.1) THEN
         iflag=0
         GO TO 510
       ENDIF
 C
 C *** REJECT LAST STEP IF ER.LT.ONE OR ISTART.EQ.1
 C
   220 CONTINUE
       IF (err.GT.one.OR.istart.EQ.1) THEN
         t=tdum
         DO 610 i=1,n
   610   x(i)=xdum(i)
         nrej=nrej+1
         nstp=nstp-1
         GO TO 100
       ENDIF
 C
 C *** CHECK FOR TIME TO CALL USEROP
 C
   510 CONTINUE
       IF (sdt.GT.zero.AND.t.GE.tpr) GO TO 200
       IF (sdt.LT.zero.AND.t.LE.tpr) GO TO 200
       GO TO 100
   200 CONTINUE
 C      CALL USEROP(T,X,N,DT,NREJ,NSTP,ISTART,IFLAG)
       IF (istart.EQ.1) GO TO 220
       GO TO 100
   900 CONTINUE
 C
 C *** CALL USEROP THE LAST TIME
 C
 C      CALL USEROP(T,X,N,DT,NREJ,NSTP,ISTART,IFLAG)
       RETURN
       END
#define f
Definition: l1_czcs.c:696
const
subroutine const(NGAUSS, NMAX, MMAX, P, X, W, AN, ANN, S, SS, NP, EPS)
Definition: ampld.lp.f:924
der
subroutine der(T, X, DX)
Definition: der.f:2
fact
subroutine fact
Definition: tmd.lp.f:1161
rk78
subroutine rk78(DER, T, TOUT, N, X, H, RELERR, ABSERR, SP)
Definition: rk78.f:2