C 
C SUBROUTINE 'RKGS' FROM THE IBM SCIENTIFIC SUBROUTINE PACKAGE.
C
C NOTE:  TO CONFORM WITH THE FORTRAN77 STANDARD, DUMMY ARRAY DIMENSIONS
C        (1) HAVE BEEN CHANGED TO (*).
C
C     ..................................................................
C
C        SUBROUTINE RKGS
C
C        PURPOSE
C           TO SOLVE A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL
C           EQUATIONS WITH GIVEN INITIAL VALUES.
C
C        USAGE
C           CALL RKGS (PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)
C           PARAMETERS FCT AND OUTP REQUIRE AN EXTERNAL STATEMENT.
C
C        DESCRIPTION OF PARAMETERS
C           PRMT   - AN INPUT AND OUTPUT VECTOR WITH DIMENSION GREATER
C                    OR EQUAL TO 5, WHICH SPECIFIES THE PARAMETERS OF
C                    THE INTERVAL AND OF ACCURACY AND WHICH SERVES FOR
C                    COMMUNICATION BETWEEN OUTPUT SUBROUTINE (FURNISHED
C                    BY THE USER) AND SUBROUTINE RKGS. EXCEPT PRMT(5)
C                    THE COMPONENTS ARE NOT DESTROYED BY SUBROUTINE
C                    RKGS AND THEY ARE
C           PRMT(1)- LOWER BOUND OF THE INTERVAL (INPUT),
C           PRMT(2)- UPPER BOUND OF THE INTERVAL (INPUT),
C           PRMT(3)- INITIAL INCREMENT OF THE INDEPENDENT VARIABLE
C                    (INPUT),
C           PRMT(4)- UPPER ERROR BOUND (INPUT). IF ABSOLUTE ERROR IS
C                    GREATER THAN PRMT(4), INCREMENT GETS HALVED.
C                    IF INCREMENT IS LESS THAN PRMT(3) AND ABSOLUTE
C                    ERROR LESS THAN PRMT(4)/50, INCREMENT GETS DOUBLED.
C                    THE USER MAY CHANGE PRMT(4) BY MEANS OF HIS
C                    OUTPUT SUBROUTINE.
C           PRMT(5)- NO INPUT PARAMETER. SUBROUTINE RKGS INITIALIZES
C                    PRMT(5)=0. IF THE USER WANTS TO TERMINATE
C                    SUBROUTINE RKGS AT ANY OUTPUT POINT, HE HAS TO
C                    CHANGE PRMT(5) TO NON-ZERO BY MEANS OF SUBROUTINE
C                    OUTP. FURTHER COMPONENTS OF VECTOR PRMT ARE
C                    FEASIBLE IF ITS DIMENSION IS DEFINED GREATER
C                    THAN 5. HOWEVER SUBROUTINE RKGS DOES NOT REQUIRE
C                    AND CHANGE THEM. NEVERTHELESS THEY MAY BE USEFUL
C                    FOR HANDING RESULT VALUES TO THE MAIN PROGRAM
C                    (CALLING RKGS) WHICH ARE OBTAINED BY SPECIAL
C                    MANIPULATIONS WITH OUTPUT DATA IN SUBROUTINE OUTP.
C           Y      - INPUT VECTOR OF INITIAL VALUES.  (DESTROYED)
C                    LATERON Y IS THE RESULTING VECTOR OF DEPENDENT
C                    VARIABLES COMPUTED AT INTERMEDIATE POINTS X.
C           DERY   - INPUT VECTOR OF ERROR WEIGHTS.  (DESTROYED)
C                    THE SUM OF ITS COMPONENTS MUST BE EQUAL TO 1.
C                    LATERON DERY IS THE VECTOR OF DERIVATIVES, WHICH
C                    BELONG TO FUNCTION VALUES Y AT A POINT X.
C           NDIM   - AN INPUT VALUE, WHICH SPECIFIES THE NUMBER OF
C                    EQUATIONS IN THE SYSTEM.
C           IHLF   - AN OUTPUT VALUE, WHICH SPECIFIES THE NUMBER OF
C                    BISECTIONS OF THE INITIAL INCREMENT. IF IHLF GETS
C                    GREATER THAN 10, SUBROUTINE RKGS RETURNS WITH
C                    ERROR MESSAGE IHLF=11 INTO MAIN PROGRAM. ERROR
C                    MESSAGE IHLF=12 OR IHLF=13 APPEARS IN CASE
C                    PRMT(3)=0 OR IN CASE SIGN(PRMT(3)).NE.SIGN(PRMT(2)-
C                    PRMT(1)) RESPECTIVELY.
C           FCT    - THE NAME OF AN EXTERNAL SUBROUTINE USED. THIS
C                    SUBROUTINE COMPUTES THE RIGHT HAND SIDES DERY OF
C                    THE SYSTEM TO GIVEN VALUES X AND Y. ITS PARAMETER
C                    LIST MUST BE X,Y,DERY. SUBROUTINE FCT SHOULD
C                    NOT DESTROY X AND Y.
C           OUTP   - THE NAME OF AN EXTERNAL OUTPUT SUBROUTINE USED.
C                    ITS PARAMETER LIST MUST BE X,Y,DERY,IHLF,NDIM,PRMT.
C                    NONE OF THESE PARAMETERS (EXCEPT, IF NECESSARY,
C                    PRMT(4),PRMT(5),...) SHOULD BE CHANGED BY
C                    SUBROUTINE OUTP. IF PRMT(5) IS CHANGED TO NON-ZERO,
C                    SUBROUTINE RKGS IS TERMINATED.
C           AUX    - AN AUXILIARY STORAGE ARRAY WITH 8 ROWS AND NDIM
C                    COLUMNS.
C
C        REMARKS
C           THE PROCEDURE TERMINATES AND RETURNS TO CALLING PROGRAM, IF
C           (1) MORE THAN 10 BISECTIONS OF THE INITIAL INCREMENT ARE
C               NECESSARY TO GET SATISFACTORY ACCURACY (ERROR MESSAGE
C               IHLF=11),
C           (2) INITIAL INCREMENT IS EQUAL TO 0 OR HAS WRONG SIGN
C               (ERROR MESSAGES IHLF=12 OR IHLF=13),
C           (3) THE WHOLE INTEGRATION INTERVAL IS WORKED THROUGH,
C           (4) SUBROUTINE OUTP HAS CHANGED PRMT(5) TO NON-ZERO.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           THE EXTERNAL SUBROUTINES FCT(X,Y,DERY) AND
C           OUTP(X,Y,DERY,IHLF,NDIM,PRMT) MUST BE FURNISHED BY THE USER.
C
C        METHOD
C           EVALUATION IS DONE BY MEANS OF FOURTH ORDER RUNGE-KUTTA
C           FORMULAE IN THE MODIFICATION DUE TO GILL. ACCURACY IS
C           TESTED COMPARING THE RESULTS OF THE PROCEDURE WITH SINGLE
C           AND DOUBLE INCREMENT.
C           SUBROUTINE RKGS AUTOMATICALLY ADJUSTS THE INCREMENT DURING
C           THE WHOLE COMPUTATION BY HALVING OR DOUBLING. IF MORE THAN
C           10 BISECTIONS OF THE INCREMENT ARE NECESSARY TO GET
C           SATISFACTORY ACCURACY, THE SUBROUTINE RETURNS WITH
C           ERROR MESSAGE IHLF=11 INTO MAIN PROGRAM.
C           TO GET FULL FLEXIBILITY IN OUTPUT, AN OUTPUT SUBROUTINE
C           MUST BE FURNISHED BY THE USER.
C           FOR REFERENCE, SEE
C           RALSTON/WILF, MATHEMATICAL METHODS FOR DIGITAL COMPUTERS,
C           WILEY, NEW YORK/LONDON, 1960, PP.110-120.
C
C     ..................................................................
C
      SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)
C
C
      DIMENSION Y(*),DERY(*),AUX(8,*),A(4),B(4),C(4),PRMT(*)
      DO 1 I=1,NDIM
    1 AUX(8,I)=.06666667*DERY(I)
      X=PRMT(1)
      XEND=PRMT(2)
      H=PRMT(3)
      PRMT(5)=0.
      CALL FCT(X,Y,DERY)
C
C     ERROR TEST
      IF(H*(XEND-X))38,37,2
C
C     PREPARATIONS FOR RUNGE-KUTTA METHOD
    2 A(1)=.5
      A(2)=.2928932
      A(3)=1.707107
      A(4)=.1666667
      B(1)=2.
      B(2)=1.
      B(3)=1.
      B(4)=2.
      C(1)=.5
      C(2)=.2928932
      C(3)=1.707107
      C(4)=.5
C
C     PREPARATIONS OF FIRST RUNGE-KUTTA STEP
      DO 3 I=1,NDIM
      AUX(1,I)=Y(I)
      AUX(2,I)=DERY(I)
      AUX(3,I)=0.
    3 AUX(6,I)=0.
      IREC=0
      H=H+H
      IHLF=-1
      ISTEP=0
      IEND=0
C
C
C     START OF A RUNGE-KUTTA STEP
    4 IF((X+H-XEND)*H)7,6,5
    5 H=XEND-X
    6 IEND=1
C
C     RECORDING OF INITIAL VALUES OF THIS STEP
    7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)
      IF(PRMT(5))40,8,40
    8 ITEST=0
    9 ISTEP=ISTEP+1
C
C
C     START OF INNERMOST RUNGE-KUTTA LOOP
      J=1
   10 AJ=A(J)
      BJ=B(J)
      CJ=C(J)
      DO 11 I=1,NDIM
      R1=H*DERY(I)
      R2=AJ*(R1-BJ*AUX(6,I))
      Y(I)=Y(I)+R2
      R2=R2+R2+R2
   11 AUX(6,I)=AUX(6,I)+R2-CJ*R1
      IF(J-4)12,15,15
   12 J=J+1
      IF(J-3)13,14,13
   13 X=X+.5*H
   14 CALL FCT(X,Y,DERY)
      GOTO 10
C     END OF INNERMOST RUNGE-KUTTA LOOP
C
C
C     TEST OF ACCURACY
   15 IF(ITEST)16,16,20
C
C     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY
   16 DO 17 I=1,NDIM
   17 AUX(4,I)=Y(I)
      ITEST=1
      ISTEP=ISTEP+ISTEP-2
   18 IHLF=IHLF+1
      X=X-H
      H=.5*H
      DO 19 I=1,NDIM
      Y(I)=AUX(1,I)
      DERY(I)=AUX(2,I)
   19 AUX(6,I)=AUX(3,I)
      GOTO 9
C
C     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE
   20 IMOD=ISTEP/2
      IF(ISTEP-IMOD-IMOD)21,23,21
   21 CALL FCT(X,Y,DERY)
      DO 22 I=1,NDIM
      AUX(5,I)=Y(I)
   22 AUX(7,I)=DERY(I)
      GOTO 9
C
C     COMPUTATION OF TEST VALUE DELT
   23 DELT=0.
      DO 24 I=1,NDIM
   24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))
      IF(DELT-PRMT(4))28,28,25
C
C     ERROR IS TOO GREAT
   25 IF(IHLF-10)26,36,36
   26 DO 27 I=1,NDIM
   27 AUX(4,I)=AUX(5,I)
      ISTEP=ISTEP+ISTEP-4
      X=X-H
      IEND=0
      GOTO 18
C
C     RESULT VALUES ARE GOOD
   28 CALL FCT(X,Y,DERY)
      DO 29 I=1,NDIM
      AUX(1,I)=Y(I)
      AUX(2,I)=DERY(I)
      AUX(3,I)=AUX(6,I)
      Y(I)=AUX(5,I)
   29 DERY(I)=AUX(7,I)
      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)
      IF(PRMT(5))40,30,40
   30 DO 31 I=1,NDIM
      Y(I)=AUX(1,I)
   31 DERY(I)=AUX(2,I)
      IREC=IHLF
      IF(IEND)32,32,39
C
C     INCREMENT GETS DOUBLED
   32 IHLF=IHLF-1
      ISTEP=ISTEP/2
      H=H+H
      IF(IHLF)4,33,33
   33 IMOD=ISTEP/2
      IF(ISTEP-IMOD-IMOD)4,34,4
   34 IF(DELT-.02*PRMT(4))35,35,4
   35 IHLF=IHLF-1
      ISTEP=ISTEP/2
      H=H+H
      GOTO 4
C
C
C     RETURNS TO CALLING PROGRAM
   36 IHLF=11
      CALL FCT(X,Y,DERY)
      GOTO 39
   37 IHLF=12
      GOTO 39
   38 IHLF=13
   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)
   40 RETURN
      END
C
C=======================================================================
C