      SUBROUTINE GERKS(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,GERROR,
     1                YP,H,F1,F2,F3,F4,F5,YG,YGP,SAVRE,SAVAE,
     2                NFE,KOP,INIT,JFLAG,KFLAG)
C
C     FEHLBERG FOURTH(FIFTH) ORDER RUNGE-KUTTA METHOD WITH
C     GLOBAL ERROR ASSESSMENT
C
C***********************************************************************
C
C     GERKS INTEGRATES A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL
C     EQUATIONS AS DESCRIBED IN THE COMMENTS FOR GERK.  THE ARRAYS
C     YP,F1,F2,F3,F4,F5,YG AND YGP (OF DIMENSION AT LEAST NEQN) AND
C     THE VARIABLES H,SAVRE,SAVAE,NFE,KOP,INIT,JFLAG,AND KFLAG ARE USED
C     INTERNALLY BY THE CODE AND APPEAR IN THE CALL LIST TO ELIMINATE
C     LOCAL RETENTION OF VARIABLES BETWEEN CALLS. ACCORDINGLY, THEY
C     SHOULD NOT BE ALTERED. ITEMS OF POSSIBLE INTEREST ARE
C         YP - DERIVATIVE OF SOLUTION VECTOR AT T
C         H  - AN APPROPRIATE STEPSIZE TO BE USED FOR THE NEXT STEP
C         NFE- COUNTER ON THE NUMBER OF DERIVATIVE FUNCTION EVALUATIONS
C
C***********************************************************************
C
      LOGICAL HFAILD,OUTPUT
C
      DIMENSION Y(NEQN),YP(NEQN),F1(NEQN),F2(NEQN),F3(NEQN),F4(NEQN),
     1          F5(NEQN),YG(NEQN),YGP(NEQN),GERROR(NEQN)
C
      EXTERNAL F
C
C***********************************************************************
C  THE COMPUTER UNIT ROUNDOFF ERROR U IS THE SMALLEST POSITIVE VALUE
C  REPRESENTABLE IN THE MACHINE SUCH THAT  1.+ U .GT. 1.
C
      U = SPMPAR(1)
C
C***********************************************************************
C
C  REMIN IS A TOLERANCE THRESHOLD WHICH IS ALSO DETERMINED BY THE
C  INTEGRATION METHOD. IN PARTICULAR, A FIFTH ORDER METHOD WILL
C  GENERALLY NOT BE CAPABLE OF DELIVERING ACCURACIES NEAR LIMITING
C  PRECISION ON COMPUTERS WITH LONG WORDLENGTHS.
C
      REMIN = 3.E-11
C
C***********************************************************************
C
C     THE EXPENSE IS CONTROLLED BY RESTRICTING THE NUMBER
C     OF FUNCTION EVALUATIONS TO BE APPROXIMATELY MAXNFE.
C     AS SET,THIS CORRESPONDS TO ABOUT 500 STEPS.
C
      MAXNFE = 9000
C
C***********************************************************************
C
C
C     CHECK INPUT PARAMETERS
C
C
      IF (NEQN .LT. 1) GO TO 10
      IF ((RELERR .LT. 0.)  .OR.  (ABSERR .LT. 0.)) GO TO 10
      MFLAG=IABS(IFLAG)
      IF ((MFLAG .GE. 1)  .AND.  (MFLAG .LE. 7)) GO TO 20
C
C     INVALID INPUT
   10 IFLAG=7
      RETURN
C
C     IS THIS THE FIRST CALL
   20 IF (MFLAG .EQ. 1) GO TO 50
C
C     CHECK CONTINUATION POSSIBILITIES
C
      IF (T .EQ. TOUT) GO TO 10
      IF (MFLAG .NE. 2) GO TO 25
C
C     IFLAG = +2 OR -2
      IF (INIT .EQ. 0) GO TO 45
      IF (KFLAG .EQ. 3) GO TO 40
      IF ((KFLAG .EQ. 4)  .AND.  (ABSERR .EQ. 0.)) GO TO 30
      IF ((KFLAG .EQ. 5)  .AND.  (RELERR .LE. SAVRE)  .AND.
     1    (ABSERR .LE. SAVAE)) GO TO 30
      GO TO 50
C
C     IFLAG = 3,4,5,6, OR 7
   25 IF (IFLAG .EQ. 3) GO TO 40
      IF ((IFLAG .EQ. 4) .AND. (ABSERR .GT. 0.)) GO TO 45
C
C     INTEGRATION CANNOT BE CONTINUED SINCE USER DID NOT RESPOND TO
C     THE INSTRUCTIONS PERTAINING TO IFLAG=4,5,6 OR 7
   30 STOP
C
C***********************************************************************
C
C     RESET FUNCTION EVALUATION COUNTER
   40 NFE=0
      IF (MFLAG .EQ. 2) GO TO 50
C
C     RESET FLAG VALUE FROM PREVIOUS CALL
   45 IFLAG=JFLAG
C
C     SAVE INPUT IFLAG AND SET CONTINUATION FLAG VALUE FOR SUBSEQUENT
C     INPUT CHECKING
   50 JFLAG=IFLAG
      KFLAG=0
C
C     SAVE RELERR AND ABSERR FOR CHECKING INPUT ON SUBSEQUENT CALLS
      SAVRE=RELERR
      SAVAE=ABSERR
C
C     RESTRICT RELATIVE ERROR TOLERANCE TO BE AT LEAST AS LARGE AS
C     32U+REMIN TO AVOID LIMITING PRECISION DIFFICULTIES ARISING FROM
C     IMPOSSIBLE ACCURACY REQUESTS
C
      RER=AMAX1(RELERR,32.*U+REMIN)
C
      U26=26.*U
C
      DT=TOUT-T
C
      IF (MFLAG .EQ. 1) GO TO 60
      IF (INIT .EQ. 0) GO TO 65
      GO TO 80
C
C
C***********************************************************************
C
C     INITIALIZATION --
C                       SET INITIALIZATION COMPLETION INDICATOR,INIT
C                       SET INDICATOR FOR TOO MANY OUTPUT POINTS,KOP
C                       EVALUATE INITIAL DERIVATIVES
C                       COPY INITIAL VALUES AND DERIVATIVES FOR THE
C                             PARALLEL SOLUTION
C                       SET COUNTER FOR FUNCTION EVALUATIONS,NFE
C                       ESTIMATE STARTING STEPSIZE
C
   60 INIT=0
      KOP=0
C
      A=T
      CALL F(A,Y,YP)
      NFE=1
      IF (T .NE. TOUT) GO TO 65
      IFLAG=2
      RETURN
C
C
   65 INIT=1
      H=ABS(DT)
      TOLN=0.
      DO 70 K=1,NEQN
        YG(K)=Y(K)
        YGP(K)=YP(K)
        TOL=RER*ABS(Y(K))+ABSERR
        IF (TOL .LE. 0.) GO TO 70
        TOLN=TOL
        YPK=ABS(YP(K))
        IF (YPK*H**5 .GT. TOL) H=(TOL/YPK)**0.2
   70   CONTINUE
      IF (TOLN .LE. 0.) H=0.
      H=AMAX1(H,U26*AMAX1(ABS(T),ABS(DT)))
C
C
C***********************************************************************
C
C     SET STEPSIZE FOR INTEGRATION IN THE DIRECTION FROM T TO TOUT
C
   80 H=SIGN(H,DT)
C
C     TEST TO SEE IF GERK IS BEING SEVERELY IMPACTED BY TOO MANY
C     OUTPUT POINTS
C
      IF (ABS(H) .GT. 2.*ABS(DT)) KOP=KOP+1
      IF (KOP .NE. 100) GO TO 85
      KOP=0
      IFLAG=6
      RETURN
C
   85 IF (ABS(DT) .GT. U26*ABS(T)) GO TO 95
C
C     IF TOO CLOSE TO OUTPUT POINT,EXTRAPOLATE AND RETURN
C
      DO 90 K=1,NEQN
        YG(K)=YG(K)+DT*YGP(K)
   90   Y(K)=Y(K)+DT*YP(K)
      A=TOUT
      CALL F(A,YG,YGP)
      CALL F(A,Y,YP)
      NFE=NFE+2
      GO TO 300
C
C     INITIALIZE OUTPUT POINT INDICATOR
C
   95 OUTPUT= .FALSE.
C
C     TO AVOID PREMATURE UNDERFLOW IN THE ERROR TOLERANCE FUNCTION,
C     SCALE THE ERROR TOLERANCES
C
      SCALE=2./RER
      AE=SCALE*ABSERR
C
C
C***********************************************************************
C***********************************************************************
C     STEP BY STEP INTEGRATION
C
  100 HFAILD= .FALSE.
C
C     SET SMALLEST ALLOWABLE STEPSIZE
C
      HMIN=U26*ABS(T)
C
C     ADJUST STEPSIZE IF NECESSARY TO HIT THE OUTPUT POINT.
C     LOOK AHEAD TWO STEPS TO AVOID DRASTIC CHANGES IN THE STEPSIZE
C     AND THUS LESSEN THE IMPACT OF OUTPUT POINTS ON THE CODE.
C
      DT=TOUT-T
      IF (ABS(DT) .GE. 2.*ABS(H)) GO TO 200
      IF (ABS(DT) .GT. ABS(H)) GO TO 150
C
C     THE NEXT SUCCESSFUL STEP WILL COMPLETE THE INTEGRATION TO THE
C     OUTPUT POINT
C
      OUTPUT= .TRUE.
      H=DT
      GO TO 200
C
  150 H=0.5*DT
C
C
C
C***********************************************************************
C     CORE INTEGRATOR FOR TAKING A SINGLE STEP
C***********************************************************************
C     THE TOLERANCES HAVE BEEN SCALED TO AVOID PREMATURE UNDERFLOW IN
C     COMPUTING THE ERROR TOLERANCE FUNCTION ET.
C     TO AVOID PROBLEMS WITH ZERO CROSSINGS,RELATIVE ERROR IS MEASURED
C     USING THE AVERAGE OF THE MAGNITUDES OF THE SOLUTION AT THE
C     BEGINNING AND END OF A STEP.
C     THE ERROR ESTIMATE FORMULA HAS BEEN GROUPED TO CONTROL LOSS OF
C     SIGNIFICANCE.
C     TO DISTINGUISH THE VARIOUS ARGUMENTS, H IS NOT PERMITTED
C     TO BECOME SMALLER THAN 26 UNITS OF ROUNDOFF IN T.
C     PRACTICAL LIMITS ON THE CHANGE IN THE STEPSIZE ARE ENFORCED TO
C     SMOOTH THE STEPSIZE SELECTION PROCESS AND TO AVOID EXCESSIVE
C     CHATTERING ON PROBLEMS HAVING DISCONTINUITIES.
C     TO PREVENT UNNECESSARY FAILURES, THE CODE USES 9/10 THE STEPSIZE
C     IT ESTIMATES WILL SUCCEED.
C     AFTER A STEP FAILURE, THE STEPSIZE IS NOT ALLOWED TO INCREASE FOR
C     THE NEXT ATTEMPTED STEP. THIS MAKES THE CODE MORE EFFICIENT ON
C     PROBLEMS HAVING DISCONTINUITIES AND MORE EFFECTIVE IN GENERAL
C     SINCE LOCAL EXTRAPOLATION IS BEING USED AND THE ERROR ESTIMATE
C     MAY BE UNRELIABLE OR UNACCEPTABLE WHEN A STEP FAILS.
C***********************************************************************
C
C
C     TEST NUMBER OF DERIVATIVE FUNCTION EVALUATIONS.
C     IF OKAY,TRY TO ADVANCE THE INTEGRATION FROM T TO T+H
C
  200 IF (NFE .LE. MAXNFE) GO TO 220
C
C     TOO MUCH WORK
      IFLAG=3
      KFLAG=3
      RETURN
C
C     ADVANCE AN APPROXIMATE SOLUTION OVER ONE STEP OF LENGTH H
C
  220 CALL FEHL(F,NEQN,YG,T,H,YGP,F1,F2,F3,F4,F5,F1)
      NFE=NFE+5
C
C     COMPUTE AND TEST ALLOWABLE TOLERANCES VERSUS LOCAL ERROR ESTIMATES
C     AND REMOVE SCALING OF TOLERANCES. NOTE THAT RELATIVE ERROR IS
C     MEASURED WITH RESPECT TO THE AVERAGE OF THE MAGNITUDES OF THE
C     SOLUTION AT THE BEGINNING AND END OF THE STEP.
C
      EEOET=0.
      DO 250 K=1,NEQN
        ET=ABS(YG(K))+ABS(F1(K))+AE
        IF (ET .GT. 0.) GO TO 240
C
C       INAPPROPRIATE ERROR TOLERANCE
        IFLAG=4
        KFLAG=4
        RETURN
C
  240   EE=ABS((-2090.*YGP(K)+(21970.*F3(K)-15048.*F4(K)))+
     1                        (22528.*F2(K)-27360.*F5(K)))
  250   EEOET=AMAX1(EEOET,EE/ET)
C
      ESTTOL=ABS(H)*EEOET*SCALE/752400.
C
      IF (ESTTOL .LE. 1.) GO TO 260
C
C
C     UNSUCCESSFUL STEP
C                       REDUCE THE STEPSIZE , TRY AGAIN
C                       THE DECREASE IS LIMITED TO A FACTOR OF 1/10
C
      HFAILD= .TRUE.
      OUTPUT= .FALSE.
      S=0.1
      IF (ESTTOL .LT. 59049.) S=0.9/ESTTOL**0.2
      H=S*H
      IF (ABS(H) .GT. HMIN) GO TO 200
C
C     REQUESTED ERROR UNATTAINABLE AT SMALLEST ALLOWABLE STEPSIZE
      IFLAG=5
      KFLAG=5
      RETURN
C
C
C     SUCCESSFUL STEP
C                        STORE ONE-STEP SOLUTION YG AT T+H
C                        AND EVALUATE DERIVATIVES THERE
C
  260 TS=T
      T=T+H
      DO 270 K=1,NEQN
  270   YG(K)=F1(K)
      A=T
      CALL F(A,YG,YGP)
      NFE=NFE+1
C
C                        NOW ADVANCE THE Y SOLUTION OVER TWO STEPS OF
C                        LENGTH H/2 AND EVALUATE DERIVATIVES THERE
C
      HH=0.5*H
      CALL FEHL(F,NEQN,Y,TS,HH,YP,F1,F2,F3,F4,F5,Y)
      TS=TS+HH
      A=TS
      CALL F(A,Y,YP)
      CALL FEHL(F,NEQN,Y,TS,HH,YP,F1,F2,F3,F4,F5,Y)
      A=T
      CALL F(A,Y,YP)
      NFE=NFE+12
C
C
C                       CHOOSE NEXT STEPSIZE
C                       THE INCREASE IS LIMITED TO A FACTOR OF 5
C                       IF STEP FAILURE HAS JUST OCCURRED, NEXT
C                          STEPSIZE IS NOT ALLOWED TO INCREASE
C
      S=5.
      IF (ESTTOL .GT. 1.889568E-4) S=0.9/ESTTOL**0.2
      IF (HFAILD) S=AMIN1(S,1.)
      H=SIGN(AMAX1(S*ABS(H),HMIN),H)
C
C***********************************************************************
C     END OF CORE INTEGRATOR
C***********************************************************************
C
C
C
C     SHOULD WE TAKE ANOTHER STEP
C
      IF (OUTPUT) GO TO 300
      IF (IFLAG .GT. 0) GO TO 100
C
C***********************************************************************
C***********************************************************************
C
C
C     INTEGRATION SUCCESSFULLY COMPLETED
C
C     ONE-STEP MODE
      IFLAG=-2
      GO TO 400
C
C     INTERVAL MODE
  300 T=TOUT
      IFLAG=2
C
  400 DO 450 K=1,NEQN
  450   GERROR(K)=(YG(K)-Y(K))/31.
C
      RETURN
      END
