text:
      subroutine SIVA(TSPECS, Y, F, KORD, NEQ, SIVAF, SIVAO, IDIMT,
     1   IDIMY, IDIMF, IDIMK, IOPT)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
c>> 2010-02-20 SIVA  Krogh  Fixed calling SIVAOP with array other than F
c>> 2009-11-03 SIVA  Krogh  Added option 11, more variables initialized.
c>> 2009-10-30 SIVA  Krogh  Gave KSSTRT and ROBND initial values.
c>> 2009-10-30 SIVA  Krogh  Fixed reference to undefined location in F.
c>> 2009-10-21 SIVA  Krogh  Got rid of NaN in diag. print when LSC=3.
c>> 2009-10-15 SIVA  Krogh  A few changes on how noise is handled.
c>> 2002-11-12 SIVA  Krogh  Fixed problem integrating to final output pt
c>> 2002-08-29 SIVA  Krogh  Added test for invalid HMIN/HMAX.
c>> 2002-07-26 SIVA  Krogh  Added KOUTKO to fully support Option 10.
c>> 2002-05-14 SIVA  Krogh  Fix starting prob. for Option 18.
c>> 2002-05-13 SIVA  Krogh  Put exponent letter in  numbers missing them
c>> 2002-05-12 SIVA  Krogh  Added error message for bad option 5 usage.
c>> 2001-09-07 SIVA  Krogh  Changes to allow user tol on G-Stops.
C>> 2001-05-25 SIVA  Krogh  Minor change for making .f90 version.
c>> 2001-05-18 SIVA  Krogh  Less computing with no error test
c>> 2001-05-17 SIVA  Krogh  Fixed so with no error test can't start dump
c>> 2001-04-24 SIVA  Krogh  Inserted comments from ivacom.
c>> 2000-12-01 SIVA  Krogh  Removed (some of) unused C1, MAXSTF, METEXT.
c>> 1999-12-28 SIVA  Krogh  Saved S in SIVACR for output consistency.
c>> 1999-08-19 SIVA  Krogh  Removed superfluous test above label 3520.
c>> 1997-04-22 SIVA  Krogh  Got rid of assigned go to's. F=0 if diag.
c>> 1996-08-26 SIVA  Krogh  Initialize F to 0 if dumping solution.
c>> 1996-08-23 SIVA  Krogh  Print TN not TSPECS(1) in error messages.
c>> 1996-05-30 SIVA  Krogh  Changed DERIVS/OUTPUT to  SIVAF/SIVAO.
c>> 1996-04-27 SIVA  Krogh  Changes to use .C. and C%%.
c>> 1996-03-30 SIVA  Krogh  Added external statement.
C>> 1996-03-25 SIVA  Krogh  Introduced TEXT1 to comply with F77.
c>> 1996-02-27 SIVA  Krogh  Fixed so DUMP not affected by ignored eqs.
c>> 1995-12-18 SIVA  Krogh  Fixed so no solution dump on 0 length integ.
c>> 1995-11-09 SIVA  Krogh  Fixed so char. data at col. 72 is not ' '.
c>> 1995-06-19 SIVA  Krogh  Fixed prob. with discon. just after restart.
c>> 1995-05-09 SIVA  Krogh  Fixed G-Stop/discontinuity code interaction
C>> 1995-04-26 SIVA  Krogh  Use KQMAXS instead of KQMAXI when LDIS>1000.
C>> 1995-04-26 SIVA  Krogh  Keep current KQL on discontinutiy.
C>> 1994-12-16 SIVA  Krogh  Fixed option 12 with K12 < 0.
C>> 1994-11-11 SIVA  Krogh  Declared all vars.
c>> 1994-11-02 SIVA  Krogh  Changes to use M77CON
c>> 1994-09-08 SIVA  Krogh  Added CHGTYP code.
c>> 1994-07-11 SIVA  Krogh  Fix to get same state with/without var. eqs.
c>> 1994-03-07 SIVA  Krogh  Allow larger order in single precision.
c>> 1994-01-14 SIVA  Krogh  Minor change to allow changing TFINAL.
c>> 1993-04-27 SIVA  Krogh  Additions for Conversion to C.
c>> 1993-04-12 SIVA  Krogh  Converted to use slightly altered MESS.
c>> 1993-04-12 SIVA  Krogh  Fixed LSC so sol. saved when HMAX is small.
c>> 1992-10-13 SIVA  Krogh  Fixed G-Stop/discontinuity code interaction.
c>> 1992-09-21 SIVA  Krogh  Fixed bug in discontinuity code.
c>> 1992-09-09 SIVA  Krogh  Fixed bug - Var. Eqs. with discontinuities.
c>> 1992-08-07 SIVA  Krogh  Storage map printed only if option 10 .ne. 0
c>> 1992-07-16 SIVA  Krogh  Restored correct discontinuity code.
c>> 1992-06-16 SIVA  Krogh  Eliminate reuse of storage for option 12.
c>> 1992-04-08 SIVA  Krogh  Removed unused labels, 1020, 2120.
c>> 1992-03-30 SIVA  Krogh  Fixed bug in SIVAOP error message.
c>> 1992-03-12 SIVA  Krogh  Simplified SIVABU, more digits in B's.
c>> 1992-01-16 SIVA  Krogh  Fixed minor bug in error messages.
c>> 1991-12-03 SIVA  Krogh  Major change for improved error checks.
c>> 1991-06-17 SIVA  Krogh  Fixed bug in checking storage allocation.
c>> 1991-04-11 SIVA  Krogh  Fixed minor bug re. option 12 in SIVAOP.
c>> 1991-03-28 SIVA  Krogh  Removed check at label 650 for KORD2I<0.
c>> 1991-02-08 SIVA  Krogh  Changed some floats to generics
c>> 1990-11-08 SIVA  Krogh  Fixed bug on TSPECS on discon.
c>> 1990-09-14 SIVA  Krogh  Fixed bug when discon. and sol. save.
c>> 1990-09-13 SIVA  Krogh  Increased dimension of BETA by 1.
c>> 1990-09-13 SIVA  Krogh  Added one more poss. on rel. error test.
c>> 1990-09-11 SIVA  Krogh  Recent change messed up getting dump output.
c>> 1990-06-05 SIVA  Krogh  Fixed bug in noise test, comments in IVACOM.
c>> 1990-05-08 SIVA  Krogh  Fixed new bug when TMARK hit in SIVAG.
c>> 1990-04-17 SIVA  Krogh  Fixed minor problem in SIVAIN error msg.
c>> 1990-04-10 SIVA  Krogh  Fixed interaction between discon. & dump.
c>> 1990-03-23 SIVA  Krogh  Fixed bug on option "-2", see 1989-12-07.
c>> 1990-03-20 SIVA  Krogh  Fixed rarely occuring loop.
c>> 1990-01-29 SIVA  Krogh  Removed unneeded labels.
c>> 1989-12-14 SIVA  Krogh  Saved common block SIVAEV.
c>> 1989-12-07 SIVA  Krogh  Added option "2" to SIVAOP.
c>> 1989-11-09 SIVA  Krogh  Made GG a save var. in SIVAHC
c>> 1989-08-21 SIVA  Krogh  Fix out of bounds ref. to V in SIVABU
c>> 1989-07-26 SIVA  Krogh  Fix bug in initial dim. check
c>> 1989-07-21 SIVA  Krogh  Code for integrating discontinuities
c>> 1987-12-07 SIVA  Krogh  Initial code.
c
c--S replaces "?": ?IVA,?IVAA,?IVABU,?IVACO,?IVACR,?IVAEV,?IVAF,?IVAHC,
c-- & ?IVAG,?IVAIN,?IVAMC,?IVAO,?IVAOP,?IVAPR,?IVASC,?IVACE,?IVAIE,
c-- & ?IVAPE,?MESS
c
c
c Note a "*" at the start of a name is used to indicate "D" for the
c double precision version and "S" for the single precision version.
c
c When converting between precisions, don't forget to change the value
c of KDIM set in parameter statements in a variety of routines, and to
c adjust comments for the data statements associated with EIBND in
c *IVACR, and B in *IVAHC.
c
c Entries
c  *IVA    Main entry for starting the package.
c  *IVAA   Main program inside the package, calls the other routines,
c          and does checks for output, and noise.  Called by the user
c          if reverse communication is used.
c  *IVABU  Back ups the solution to the current base time, if a step
c          that has been started must be taken over for some reason.
c  *IVACO  Called by user to get certain information from the common
c          blocks.
c  *IVACR  Corrects the solution, estimates errors, and selects order.
c  *IVADB  Subroutine to assist in debugging codes.  Called by user to
c          get a formatted list of all the variables used in the
c          integration.  Not required in usual case.
c  *IVADE  Needed only for delay differential equations.  This is called
c          by the user from the derivative subprogram.
c  *IVAG   Required only if the user has G-Stops, i.e. places to call
c          his output subroutine where certain functions have zeroes.
c  *IVAHC  Compute coefficients that depend on the step size history.
c  *IVAIN  Used to interpolate to arbitrary points.
c  *IVAOP  Used to process user option requests.
c  *IVAPR  Used to update the differences and to predict the solution
c          at the end of the current step.
c
c External Routines
c  *1MACH  Not used in the Fortran 95 version.  ("*" is "D" for double
c          and "R" for single precision.) This returns constants that
c          depend on the floating point arithmetic.  Input arguments of
c          1 to 4 give respectively:  underflow limit, overflow limit,
c          smallest relative difference between two floating point
c          numbers, and the largest relative difference between two
c          floating point numbers.
c DERIVS (formal) Name of subroutine to be called for computing

c  OPTCHK  Used in checking storage allocation.
c  *MESS   Used to output error messages and diaganostic messages.
c          (Just MESS if no floating point is output.)
c  *ZERO   Called only if *IVAG is used.  Iterates to find zeros of
c          arbitrary (continuous) functions.
c
c Common blocks -- As a left over from the distant past, some variables
c   are in common so that they would be saved.
c  *IVAEV  Holds variables that depend on the environment.
c  *IVAMC  The main common block for the package.
c  *IVASC  The secondary common block for the package.  This contains
c          variables that are required for doing interpolation and is
c          separate to simplify saving the variables that are required
c          when the solution is being dumped (saved).
c
c Common variables and local variables
c ALPHA  (*IVAMC) Array with I-th entry = (current step size) / XI(I).
c   Used in computing integration coefficients.
c B      (*IVAHC) Array used to get started on computing integration
c   coefficients.  B(K) = 1. / (K*(K+1))
c BAKMIN (*IVADE) The largest delay at the initial point.
c BETA   (*IVAMC) Array with I-th entry = product (K=1,I-1) of
c   (current (XI(K)) / XI(K) from previous step),  BETA(1)=1.  Used in
c    updating the difference tables.
c C      (*IVAIN) Array used to hold integration/interpolation coeffs.
c C0     Parameter = 0. (in *IVAA,DE,CR,A,G,HC,IN,OP,PR)
c C1     Parameter = 1. (in *IVA,A,CR,DA,HC,IN,OP)
c C10    Parameter = 10. (in *IVAA,CR,OP)
c C1000  Parameter = 1000. (in *IVACR)
c C16    Parameter = 16. (in *IVAA,OP)
c C1M3   Parameter = .001 (in *IVAA)
c C1M5   Parameter = .00001 (in *IVAA)
c C1P125 Parameter = 1.125 (in *IVAA,HC,OP)
c C1P3   Parameter = 1.3 (in *IVAA)
c C1P4   Parameter = 1.4 (in *IVACR)
c C2     Parameter = 2. (in *IVAA,DE,BU,CR,IN,OP)
c C20    Parameter = 20. (in *IVACR)
c C2P5M3 Parameter = .0025 (in *IVAA)
c C4     Parameter = 4. (in *IVACR,OP)
c C40    Parameter = 40. (in *IVACR)
c C4096  Parameter = 4096. (in *IVAA)
c C6     Parameter = 6. (in *IVAA)
c C8M3   Parameter = .008 (in *IVAA)
c CM2    Parameter = -2. (in *IVACR)
c CM8    Parameter = -8. (in *IVACR)
c CMP5   Parameter = -.5 (in *IVACR)
c CMP75  Parameter = -.75 (in *IVAOP)
c CP0625 Parameter = .0625 (in *IVAA)
c CP1    Parameter = .1 (in *IVAA,CR,DA,HC)
c CP125  Parameter = .125 (in *IVACR)
c CP25   Parameter = .25 (in *IVAA,CR,DE,OP)
c CP3    Parameter = .3 (in *IVAA,OP)
c CP4    Parameter = .4 (in *IVAA)
c CP5    Parameter = .5 (in *IVAA,CR,DA,DE,HC,OP)
c CP5625 Parameter = .5625 (in *IVAHC)
c CP625  Parameter = .625 (in *IVAOP)
c CP75   Parameter = .75 (in *IVACR,OP)
c CP8    Parameter = .8 (in *IVACR)
c CP875  Parameter = .875 (in *IVAA, OP)
c CP9    Parameter = .9 (in *IVAOP)
c CP9375 Parameter = .9375 (in *IVACR)
c CQ3125 Parameter = .03125 (in *IVACR)
c CRBQI  Parameter = .421875 (in *IVAHC)  Initial val for computing RBQ.
c CSUM   (*IVAIN) Array used to contain partial sums of the integration
c   coefficients.  This is used to corrrect for a difference table that
c   has not yet been updated.
c D      (*IVAMC) Array to be used later to store coefficients for
c   integrating stiff equations.
c   derivatives.  Not used if option 13 is set.
c DISADJ (*IVAA) Value of stepsize when discontinuity is indicated.
c DNOISE (*IVAMC) Used in determining if noise is limiting the
c   precision.  It is usually |highest difference used in correcting|
c   of the equation with the largest error estimate.
c DS     (*IVAMC) Array to be used later to store coefficients for
c   estimating errors when integrating stiff equations.
c DVC2   (*IVADB) Array used for output of variables HC to TOUT in
c   common block *IVAMC.
c E      (*IVACR) (Estimated error) / (Requested accuracy)
c EAVE   (*IVAMC) This is a weighted average of past values of EIMAX.
c   It is adjusted to account for expected changes due to step changes.
c EEPS10 (*IVAEV) = 10. * (machine epsilon).
c EEPS16 (*IVAEV) = 16. * (machine epsilon).
c EEPS2  (*IVAEV) =  2. * (machine epsilon).
c EEPT75 (*IVAEV) = (machine epsilon) ** (.75)
c EI     (*IVACR) Estimate for what E would be if step size increased.
c EIBND  (*IVACR) Array containing limits on the estimated error with
c   the stepsize increased.  This array tends to make the code a little
c   more conservative on step size increases at low order.
c EIMAX  (*IVAMC) Estimate of (error estimate / error requested) if the
c   step size should be increased.
c EIMIN  (*IVAMC) An error estimate is small enough to allow a step
c   increase if the estimate of ((error with the step size increased) /
c   (error requested)) is less than EIMIN.
c EIMINO (*IVAA) Set to C8M3 and never changed.  When step size is being
c   reduced if EIMIN .le. EIMINO then the reduction factor is set to
c   CP875.  This variable could be a parameter.
c EMAX   (*IVAMC) Largest value computed for (error estimate) / (error
c   requested).
c EOVEP2 (*IVAEV) = EEPS2 * (largest floating point number).
c EPS    (*IVACR) Current absolute error tolerance.  Also used for
c   temporary storage when computing the desired value of EPS.
c ERCOEF (*IVACR) (Error coefficient from formula) / EPS
c EREP   (*IVAMC) If EMAX > EREP, a step is repeated.  Ordinarily
c   this has the value .3.  This is set < 0 if the error tolerance is
c   specified improperly, and is set to a large value if the user
c   requests complete control over the step size.  EREP is also set
c   < 0 after a user specified discontinuity.
c EROV10 (*IVAEV) = 10. / (largest floating point number).
c ETA    (*IVAIN) Array used in computing integration/interp. coeffs.
c EVC    (*IVADB) Array used for output of variables EEPS2 to EROV10 in
c   common block *IVAEV.
c EXR    (*IVAA) Set to CP1 and never changed.  If it is estimated the
c   the (error estimate) / (error requested) on the next step will be
c   .ge. EXR then the step size is reduced.  Could be a parameter.
c F      (formal) Array used to store derivative values, the difference
c   tables, error tolerance requests, and values used by some other
c   options. (in *IVA,A,BU,CR,DA,DB,G,IN,PR)
c FDAT  (*IVAMC) Used to store data for error messages.  (Local array in
c   *IVAIN.)
c FOPT  (formal) in *IVAOP.  Passed as place to save floating point data
c   for options.  This package passes F in for FOPT when calling *IVAOP.
c G      (*IVAMC) Integration coefficients used for predicting solution.
c   G(I, J) gives the I-th coefficient for integrating a J-th order
c   differential equation.  G(1, 1) is equal to the step size.
c GAMMA  (*IVAIN) Array used in computing integration/interp. coeffs.
c GG     (*IVAHC) Array of length = max. differential equation order
c   allowed by code - 1.  GG(K) = (HH**(K+1)) / K!
c GNEW   (formal) in *IVAG.  Current value for vector function g, whose
c   zeroes are to be found.
c GOINT  (*IVACR) Used for assigned go to used in computing integration
c   coefficients.
c GOLD   (*IVAG) Previous value for element of G whose zero search is
c   active.
c GS     (*IVAMC) Integration coefficients used in estimating errors.
c GT     (formal) in *IVAG.  Previous value of GNEW.
c HC     (*IVAMC) Ratio of (new step size) / (old step size)
c HDEC   (*IVAMC) Default value to use for HC when reducing the step
c   size.  (Values closer to 1 may be used some of the time.)
c HH     Equivalenced to G(1,1) = current step size in *IVAA,CR,DA,G,HC.
c HI     (*IVAIN) Step length from the base value of the independent
c   variable for the interpolation.
c HINC   (*IVAMC) Default value to use for HC when increasing the step
c   size.  (Values closer to 1 may be used some of the time.)
c HINCC  (*IVAMC) Actual value used for default value of HC when
c   increasing the step size.  Set to HINC after start is considered
c   complete.  During the start HINCC is set to 1.125.
c HMAX   (*IVAMC) Largest value allowed for abs(step size).  Default
c   value is a very large number.
c HMAXP9 (*IVAMC) .9 * HMAX.
c HMIN   (*IVAMC) Smallest value allowed for abs(step size).  Default
c   value is 0.
c HNEW   (*IVADE) Value of step size when iterating at initial point
c   for delay differential equations.
c I      Used for temporary storage. (*IVAA,BU,CR,DA,DE,G,IN,OP,PR)
c IA     (*IVAOP) absolute value of first integer stored for an option.
c ICF    (*IVAMC) Final index for current loop in *IVACR.  Required by
c   option 18.
c ICI    (*IVAIN) Temporary index, = 0 for interpolation, 1 or 0 for
c   differentiation, and d-1, d-2, ... 0 for integration, where d is the
c   order of the differential equation.  Index of first location
c   in C() used is ICI + an offset.
c ICS    (*IVAMC) Starting index for current loop in *IVACR.
c ID     (formal) Array use to contain integer data from common.  Values
c   are returned in locations 1 to 5 as follows.
c   1    KEMAX = Index of equation with largest error estimate
c   2    KSTEP = Current step number
c   3    NUMDT = Number of differences used for each equation
c   4            Reserved for future use
c   5            Reserved for future use
c IDAT   (*IVAMC) Used to store integer for error messages.  (Also used
c   in *IVAA for temporary storage of KORD(2).  (Local array in *IVAIN.)
c IDE    (*IVADE - formal) Array used to contain past information so
c   that delays can stretch back indefinitely.  If the first location is
c   0, then any interpolations requested must be in the range of the
c   current difference tables.  At present, only the value 0 is allowed
c   in IDE(1).  This array is intended for the support of saving long
c   past histories.  IDE(2) must contain the declared dimension of WDE.
c IDEF   (*IVADE -- formal) Flag giving indicaion of what is going on.
c   = 0  User should compute derivatives and return to the main
c        integrator.
c   = 1  Code is computing additional values in order to get past data
c        necessary for starting.  User should compute derivatives and
c        call *IVADE.
c   < 0  Indicates an error condition.  If *IVADE is called without
c        changing the value of IDEF, the integration is stopped and an
c        error message printed.  Possible error flags are:
c    -1  Difference tables do not span back far enough to compute the
c        past values of Y needed.
c    -2  There is not enough space in WDE to get the required starting
c        values.
c IDIMF  (formal) Declared dimension of F().
c IDIMK  (formal) Declared dimension of KORD().
c IDIMT  (formal) Declared dimension of TSPECS().
c IDIMY  (formal) Declared dimension of Y().
c IDT    (*IVAIN) Used as a base index into the difference table.
c IFLAG  (formal in *IVAG) Used for communication with user.
c   = 1  Continue as if *IVAG was not called.
c   = 2  Check KORD(1) as one would do at start of OUTPUT if no G-Stops
c        were present. (Exit if in DERIVS.)
c   = 3  Return to the integrator.
c   = 4  Compute G and return to *IVAG.
c   = 5  A G-Stop has been found, and NSTOP gives its index.  (If NSTOP
c        < 0, the stop was an extrapolating stop.)
c   = 6  Same as 5, but requested accuracy was not met.
c   = 7  Same as 5, but there is a probable error in computing G.
c   = 8  Fatal error of some type.  (An error message has been printed.)
c IG     (*IVAG)  IG = KORD(2) on the initial entry (0 for extrapolating
c   G-Stops, and 1 for interpolating).
c IGFLG  (*IVAMC) Used primarily in *ivag, but also used in *iva to keep
c   track of the state of GSTOP calculations.
c   = -2 Extrapolatory G's initialized, but not the interpolatory.
c   = -1 Interpolatory G's initialized, but not the extrapolatory.
c   =  0 Set when integration is started or restarted, or option setting
c        GSTOP is set.
c   =  1 Iterating to find a GSTOP.
c   =  2 User told that a GSTOP was found.
c   =  3 Checking G's at point where a GSTOP was located.
c   =  4 Checking G's at a T output point.
c   =  5 Usual case, no sign change detected.
c IGSTOP (*IVAMC) IGSTOP(k) is set in *ivag to the index of the last G
c   with a 0, where k is one for an interpolatory G-Stop, and k is two
c   for an extrapolatory G-Stop.
c IGTYPE (*IVAMC) Array with two elements as for IGSTOP, but this saves
c   a flag giving the nature of convergence to the stop.
c   = 0  All known G-stops completely processed.
c   = 4  Need to compute next value while iterating.
c   = 5  Got good convergence.
c   = 6  Got convergence, but not to desired accuracy.
c   = 7  Problem in getting convergence.
c   = 8  A fatal error of some type.
c IHI    (*IVA) Last location used by the current option.
c ILGREP (*IVAMC) Used when correction to keep track of equations that
c   are to use a certain error tolerance.
c ILGROR (*IVACR) Index of last equation in the current group of
c   equations grouped for selecting integration order.
c ILOW   (*IVA) First location used by the current option.
c INCOM  (*IVADE) Array equivalenced to LDT in the common block *IVASC.
c   Used to simplify saving information in the common block.
c INCOP  (*IVAOP) Array containing data giving the amount of space in
c   IOPT used for each of the options.
c INGS   Current index for G-stop being examined in SIVAG.
c INICAS (*IVADE) Used to track the initialization for a delay equation.
c   = 1  Very beginning.
c   = 2  Getting derivative at the very beginning.
c   = 3  Getting derivatives at points prior to the initial point.
c   = 4  Getting derivative at initial point after iteration is started.
c INTCHK (*IVA) Array passed to OPTCHK containing information on storage
c   allocation.  See comments in OPTCHK for details.
c INTEG  (*IVAIN) Number of integrations being done. (<0 for
c   differentiations and =0 for interpolation.)  Also used as counter
c   when computing integration coefficients.
c        (*IVAPR) Number of integrations being done.
c INTEGS (*IVAPR) = -1 for equations that are not stiff, 0 for those
c   that are stiff.
c INTEGZ (*IVAIN) min(INTEG, 0)
c INTERP (*IVAIN) added to the usual integration order to get the order
c   to be used when interpolating: 3-KQMAXI, if HI=0; 1, if
c   |HI| > |XI(1)| and HI * XI(1) < 0; 0, otherwise -- the usual case.
c IOP10  (*IVAMC) Number of times diagnostic output is to be given when
c   leaving *ivacr (the corrector).
c IOP11  (*IVAMC) Gives current step number of the method.  Tells how
c   many of certain coefficients must be computed. (Has nothing to do
c   with options.) = min(max integ order + 1, KDIM).  Also set when
c   starting to flag that certain memory locations must be set to 0.
c IOP12  (*IVAMC) Points to location in F() where user supplied values
c   of HINC, HDEC, HMIN, and HMAX.  (0 if option 12 not used.)
c IOP13  (*IVAMC) If not zero, reverse communication will be used for
c   getting the values of derivatives.  Associated with option 13.
c IOP14  (*IVAMC) If not zero, reverse communication will be used in
c   place of calls to the output routine.  Associated with option 14.
c IOP15  (*IVAMC) If not zero, a return will be made to the user after
c   the initialization.  Associated with option 15.  This might be used
c   to overlay *iva, some of the user's code, and perhaps *ivaop.
c IOP16  (*IVAMC) Points to location in KORD() where information for
c   specifying the error tolerance is specified.  See option 16.
c IOP17  (*IVAMC) Used in initialization for option 17, afterwards this
c   cell is used by KEXIT which is equivalenced to IOP17.
c IOP18  (*IVAMC) Points to location in KORD() where information for
c   specifying a grouping of equations for derivative evaluation is
c   stored.  See option 18.
c IOP19  (*IVAMC) Points to location in KORD() where information for
c   specifying a grouping of equations for integration order control
c   is stored.  See option 19.
c IOP20  (*IVAMC) Used for option 20, gives first location in F where
c   estimated errors are to be stored.  Expected to be useful in a
c   program for solving boundary value problems using multiple shooting.
c IOP21  (*IVAMC) Was used for stiff equations option (never completely
c   coded).  The optional code still uses this (don't activate it!).
c   Now used to flag the location if F where the user has stored the
c    tolerance to use in finding G-Stops.
c IOP21S (*IVAMC) Was used for stiff equations see above.
c IOP22  (*IVAMC) Set aside for possible option for stiff equations.
c IOP3   (*IVAMC) Value set by option 3.
c   =  0 Interpolate to final point. (The default)
c   =  1 Integrate to final point.
c   = -1 Extrapolate to final point.
c IOP4   (*IVAMC) Value set by option 4.  The output routine is called
c   with KORD(1) = 4, every IOP4 steps.  (Default value for IOP4 is a
c   very large number.
c IOP5   (*IVAMC) Value provided by option 5, used to specify extra
c   output points.
c IOP6   (*IVAMC) Value provided by option 6.  If nonzero, the output
c   routine is called at the end of every step.  If > 0, there are
c   IOP6 interpolating G-Stops.
c IOP7   (*IVAMC) Value provided by option 7.  If > 0, there are K7
c   extrapolating G-Stops.
c IOP8   (*IVAMC) Value provided by option 8.  If nonzero, the output
c   routine is called with KORD(1)=8 whenever the step size is changed.
c IOP9   (*IVAMC) Value provided by option 9.  Used to specify that the
c   user wishes to save the solution.
c IOPIVA (*IVA) Used to save length of IOPT vector for error messages.
c IOPST  (*IVASC) Intended for possible use in stiff equations.
c IOPT   (formal *IVA and IVAOP) Used to specify options.
c IOPTC  (*IVAOP) In *IVAOP equivalenced so that IOPTC(3) is equivalent
c   to IOP3.
c IOPTS  (*IVAOP) Array containing the current default values to be
c   stored into IOPTC.
c IORD   (*IVACR) Index of first equation in the current group of
c   equations grouped for selecting integration order.
c IOUTKO (*IVADC) Used in *IVADI to point to KORD to keep track of
c   equation grouping for diagnostic output.
c ISVCOM (*IVADE) Used to save info. in the common block *IVASC.
c ITERS  (*IVADE) Counts iterations in starting delay differential
c   equations.  Max. value for this is arbitrarily 100.
c ITOLEP (*IVAMC) Used for temporary storage, and for the index of a
c   tolerance relative to the start of tolerances.
c IVC1   (*IVADB) Array used for output of variables IOPST to NUMDT in
c   common block *IVASC.
c IVC2   (*IVADB) Array used for output of variables ICF to NY in
c   common block *IVAMC.
c IWB    (*IVADE) Current base index for saving F values in WDE when
c   starting delay differential equations.
c IY     (*IVAMC) Used for the current index to the Y() array.  (Local
c   variable in *IVAIN used in computing IYI.)  Equivalenced to
c   IZFLAG in *IVAG.
c IYI    (*IVAIN) Y(IYI) is currently being computed.
c IYN    (*IVAIN) Y(IYN) is base Y() corresponding to Y(IYI).
c IYNI   (*IVAIN) Used as base index for computing IYN as IY is for INI.
c IYO    (*IVADE) Points to first base value of Y for current
c   interpolation when getting values for a delay differential equation.
c IZFLAG (*IVAG)  Equivalenced to IY.  Set to 0 initially, and later
c   set to the value returned by *ZERO.
c    = 0  Value set on entry at start of search.
c    = 1  Compute next g again.
c    = 2  Normal terminiation.
c    = 3  Normal termination -- error criterion not satisfied.
c    = 4  Apparent discontinuity -- no zero found.
c    = 5  Couldn't find a sign change.
c    = 6  *ZERO was called with a bad value in IZFLAG.
c J      For temporary storage. (In *IVA,A,BU,CR,DA,DB,DE,HC,IN,OP,PR)
c J1     (*IVAA & DA) Used for temporary storage.
c J2     (*IVAA) Used for temporary storage.
c JL     (*IVA) Used for checking storage.
c JLGREP (*IVACR) Contents of first location of KORD (called LGROUP in
c   *IVACR) for the current error tolerance rule.
c JLGROR (*IVACR) Contents of first location of KORD for the current
c   integration order control.
c JLIM   (*IVA) Used for checking second item in KORD list for options
c   16 and 19.
c K      For temporary storage.  (In *IVA,A,BU,CR,DA,DB,DE,HC,IN,OP,PR)
c KDIM   Parameter giving the largest number of differences supported.
c        Used in all the routines.
c KEMAX  (*IVAMC) Index associated with equation giving the largest
c   value for (estimated error) / (requested error).
c KEXIT  (*IVAMC) Equivalenced to IOP17 which is not used after
c   initialization.  Defines actions when KORD2I = -7.  (Referenced in
c   (*IVAA,DA,G).)
c   =  1  Take the step over with reduced H.
c   =  2  Take the step over.
c   =  3  Do the end of step call to OUTPUT.
c   =  4  Reset TMARK, then do same as for KEXIT = 2.
c   =  5  Reset TMARK, then do same as for KEXIT = 3.
c   =  6  Give the fatal error diagnostic.
c KFERR  (*IVA)  Temporary storage in checking for option 16.
c KGO    (*IVA)  Used to tell from whence a check is being done or an
c   error message is being written.
c   = 1 Checking an equation group for variational equations.
c   = 2 Checking an equation group for diagnostic print.
c   = 3 Checking an equation group for integration order control.
c   = 4 Checking an equation group for error control.
c   = 5 Checking an equation group for specifying ODE orders.
c   = 6 Found a problem with output type for printing.
c   = 7 Found a problem with an output group for printing.
c   = 8 Found a problem with input NEQ.
c   = 9 Order specified for the ODE's in the system is out of range.
c   =10 Option 16 was not used (an error).
c   =11 Error tolerance of 0 specified without proper flags.
c KIS    (*IVAMC) Used to check if it is time to dump the solution.
c   The check involves incrementing KIS at the end of the step, and
c   dumping the solution if KIS is 0.
c   = -1  Set in *ivacr when it is time to dump solution
c   =  0  When starting
c   =  2  After being dumped.
c   This is set to 1000 just after a user specified discontinuity, and
c   counted up from that point.
c KMARK  (*IVAMC) Identifies the type of output associated with the next
c   output point specified by TSPECS.
c KONV   (*IVADE) Counts iterations.  Test for convergence if KONV > 1.
c KORD   (formal in *IVA,A,BU,CR,DA,DB,DE,G,IN,PR) KORD(1) is used to
c   return flags for the user to test, and KORD(2) tells what routine
c   the flag is associated with.  See KORD1I and KORD2I below and the
c   write up for the program.  KORD(3) is used for communicating extra
c   information to the user in some cases.  KORD(4) to KORD(NTE+3) are
c   used for integration order for the equations, and the rest of KORD()
c   is available for user options.
c KORD1I (*IVAMC) Helps in defining the state of the integrator.
c   Frequently has the same value as KORD(1).  Meaning depends on the
c   value of KORD(2), or the value about to be assigned to KORD(2).
c   <  0  Happens when preparing to give output with extrapolation.
c   =  0  Happens when checking F at points for noise test.
c   =  1  (KORD(2)=-1)  End of integration has been reached.
c   =  1  (KORD(2)= 0)  Computing first predicted derivative.
c   =  1  (KORD(2)= 1)  Output for initial point.
c   =  2  (KORD(2)=-1)  Giving diagnostic for noise limiting precision.
c   =  2  (KORD(2)= 0)  Computing corrected derivative.
c   =  2  (KORD(2)= 1)  Output for TSPECS(3).
c   =  3  (KORD(2)=-1)  Diagnostic for step size reduction too fast.
c   =  3  (KORD(2)= 0)  Computing variational derivative.
c   =  3  (KORD(2)= 1)  Output for TSPECS(4).
c   =  4  (KORD(2)=-1)  Error, discontinuity.
c   =  4  (KORD(2)= 1)  Output for certain number of steps.
c   =  5  (KORD(2)= 0)  Get initial derivatives for stiff equations.
c   =  5  (KORD(2)= 1)  Extra output from TSPECS.
c   =  6  (KORD(2)= 1)  End of step output.
c   =  7  (KORD(2)= 0)  Evaluate G before extrapolated output point.
c   =  7  (KORD(2)= 1)  Evaluate G before extrapolated output point.
c                       (Also used when checking for other G's after
c                        finding one.)
c   =  8  (KORD(2)= 1)  Tell user step size has changed.
c   =  9  (KORD(2)= 1)  Request for user to save solution.
c   = 11  (KORD(2)=-1)  Error, step size too small at end of start.
c   = 12  (KORD(2)=-1)  Error, step size is too small.
c   = 13  (KORD(2)=-1)  Error, output points specified badly.
c   = 21  (KORD(2)=-1)  H too small to give reasonable change when added
c                       to T.
c   = 22  (KORD(2)=-1)  Error, bad tolerance.
c   = 23  (KORD(2)=-1)  Set after message for a fatal error.
c   = 24  Set on error message in *iva, along with KORD2I = -4.
c   Also used as an index into MLOC in *IVAA when an error is being
c   processsed, see MLOC below.
c KORD2I (*IVAMC) Helps in defining the state of the integrator.
c   Frequently has the same value as KORD(2).
c   = -3  Set in *ivag, to get a derivative evaluation.
c   = -2  Set in *ivag, to get another entry to OUTPUT.
c   = -1  Return to calling program, done, interrupt, or got an error.
c   =  1  Calling OUTPUT or returning to user for OUTPUT type action.
c   =  0  Calling DERIVS or returning to user for DERIVS type action.
c   = -4  Error message in *iva and in *ivaop, along with KORD1I = 24.
c   = -5  Starting
c   = -6  Starting, getting the initial derivative value or derivatives
c         for the noise test.
c   = -7  Done some extrapolation, KEXIT defines the action to take.
c         Set in *ivag to activate KEXIT action in *iva.
c   = -8  Set when user has requested adjustment of the difference
c         tables for a discontinutiy.
c KORDI  (*IVASC) Order of differential equation being integrated.  If
c   all orders are the same, this set once at the beginning.
c KOUNT   (*IVADE) Count of number of points back from the initial point
c   when solving a delay differential equation.
c KOUNTM  (*IVADE) Largest value currrently allowed for KOUNT.
c KOUNTX  (*IVADE) Largest value allowed for KOUNTM.
c KOUTKO  Used in SIVACR to track where output is wanted.
c KPRED  (*IVAMC) Value assigned to KORD1I when getting a predicted
c   derivative.  (1 used now, 5 planned for use with stiff equations.)
c KQD    (*IVACR) = max(2, integration order)
c KQDCON (*IVAMC) Number of coefficients computed with constant step
c   size for stiff equations.
c KQICON (*IVAMC) Number of coefficients computed with constant step
c   size for nonstiff equations.
c KQL    (*IVACR) Integration order at start of (*IVACR)
c KQLORD (*IVACR) Saved value of KQL when equations are grouped for
c   controlling the integration order.
c KQMAXD (*IVASC) Maximum integration order used for stiff equations.
c KQMAXI (*IVASC) Maximum integration order used for nonstiff equations.
c KQMAXS (*IVAMC) Maximum integration order for equations that have
c   some limit on the error that can be committed.
c KQMXDS (*IVAMC) Used to save KQMAXD in case step is repeated and the
c   solution must be dumped.
c KQMXI  (*IVAIN) Maximum integration order used for integration or
c   interpolation, = KQMAXI+INTERP-1.
c KQMXS  (*IVAIN) Maximum step number, = max(KQMXI, KQMAXD).
c KQMXIL (*IVAMC) Value of KQMAXI the last time integration coefficients
c   were computed.
c KQMXIP (*IVAMC) = KQMAXI + MAXINT, for computing integration coeffs.
c KQMXIS (*IVAMC) Used to save KQMAXI in case step is repeated and the
c   solution must be dumped.
c KQN    (*IVACR) Value of integration order at end of *IVACR.
c KQQ    Used for the integration order for current equation.  (Values
c   < 0 are intended for stiff equations.)  (In *IVA,BU,DA,IN,PR)
c KSC    (*IVAMC) Number of steps that have been taken with a constant
c   step size.
c KSOUT  (*IVAMC) When KSTEP reaches this value, the output routine is
c   called with KORD(1) = 4.  The default value is a very large number.
c KSSTRT (*IVAMC) Set when ending one derivative per step to KSTEP + 2.
c   Checked later in *IVAHC to decide whether to set the step changing
c   factors to their nominal values.
c KSTEP  (*IVAMC) Number of steps taken since the start of integration.
c L      Used for temporary storage.  In *IVAIN, L is the initial value
c   of LDT, except L=1 if LDT=-1, and MAXINT .ge. 0.  (Used in *IVAA,BU
c   CR,DA,DB,IN,PR.)
c LAHAG  (*IVADB) Used to get proper offset into an diagnostic message.
c LAIAG  (*IVADB) Used to get proper offset into an diagnostic message.
c LDIS   (*IVAA) Count of steps since user flagged a discontinuity.
c LDT    (*IVASC) Used to keep track of state of difference table.
c   = -5  Used only on first step to indicate that an extra iteration
c         is desired to get a firm estimate on the error.
c   = -4  Set on initialization before there is any difference table.
c   = -3  Set just after predicting, interpolation is not allowed when
c         this value is set.
c   = -2  Set when difference table is to be updated to the end of the
c         current step, but no interpolation is to be done.  (For
c         dumping the solution.)
c   =  0  Calculations for current step are complete, it is not o.k. to
c         update the difference table.
c   =  1  Difference table has been updated to the end of the current
c         step, e.g. by doing an interpolation.
c   =  2  Set when doing a special interpolation during computation of
c         derivatives.  (For delay equations.)
c LEX    (*IVAMC) Indicates how to get values at next output point:
c   = -1  Extrapolate
c   =  0  Interpolate (The usual case.)
c   =  1  Integrate to the output point, integration is not continued.
c LGO    (*IVAIN) Used as an an assigned go to.  Result is to add in
c   extra correction term when LDT has been set to 2.
c LGROUP (formal) This is a part of KORD passed into *IVACR.  The first
c   location is the start of the information on the grouping of
c   equations for error control.
c LINC   (*IVAMC) Used to indicate state of step size selection.
c   = -10 After computed derivatives at base time, after computing other
c         extra derivatives for the noise test.
c   = -9  After computed second extra derivative for noise test.
c   = -8  After computed first extra derivative for noise test.
c   = -7  Dumping the solution and then doing a user initiated restart,
c         or getting ready to compute extra derivatives for the noise
c         test.
c   = -6  Dumping the solution before a restart.
c   = -5  Set on the first step, and also set when dumping the solution
c         after a discontinuity.
c   = -4  Repeat step with no change in the step size.
c   = -3  Set when the error tolerance is set improperly.
c   = -2  User has complete control of selecting the step size.
c   = -1  Step is being repeated.
c   =  0  Step size is not to be increased on this step.
c   = k>0 Step size can be increased by HINCC**k.
c LINCD  (*IVAMC) Value of smallest k for which HINCC**k .ge. 2.
c   (=-2 if user is specifying all step size changes.)
c LINCQ  (*IVAMC) Value of smallest k for which HINCC**k .ge. 4.
c LIOPT  (*IVAOP) Value of the last index in IOPT on the last call.
c   Used so *IVA can print IOPT in error messages.
c LL     (*IVACR) Temporary variable used when equations are grouped
c   for integration order control.
c LNOTM1 (*IVAIN) Logical variable = L .ne. -1.  If LNOTM1 is true,
c   storage in Y() is different in some way lost to antiquity.  Such
c   a case can only arise in the case of stiff equations.
c LOCF1  (*IVADB) Gives packed data needed for output of tables by the
c   message processor MESS.  See comments there under METABL for defs.
c LOCF2  (*IVADB) As for LOCF1 above.
c LOCM   (*IVAA) Parameter = 32*256, used to unpack integers stored
c   in MLOC for use in error message processing.
c LPRINT (formal, *IVADB) Defines how much printing is to be done in
c   *IVADB.  Let |LPRINT| = 10*N1 + N2     (N1,N2 digits)
c    N1=1   Do not print any variables external to the integrator.
c    N1=2   Print  tspecs, current y, past y, current f, all pertinent
c           contents of KORD, and TOL.
c    N1=3   Above + difference tables up to highest difference used.
c    N1=4   Same as N1=1 + all in storage allocated for differences.
c    N2=1   Do not print any variables internal to the integrator.
c    N2=2   Print all scalar variables in interpolation common block.
c    N2=3   Above + all scalar variables in main integ. common block.
c    N2=4   Same as N1=3 + all used in arrays XI,BETA,ALPHA, first
c           column of G, GS,RBQ,SIGMA
c    N2=5   Same as N1=4 + all used in arrays G,D,DS,V
c LSC    (*IVAMC) Indicates if starting or if noise may be present.
c   =k<0 -k steps have been taken for which noise appears to be limiting
c        the precision.
c   = 0  Usual case
c   = 1  Doing 1 derivative per step after initial part of start.
c   = 2  Used as flag that it is time to set LSC=0.
c   = 3  Third step, hold the order constant.
c   = 4  Second step, increase orders from 2 to 3.
c   = 5  First step, third time through the first step (if required).
c   = 6  First step, second time through.
c   = 7  First step, first time through.
c   = 8  Set on initialization.
c LTXT?? Names of this form are used in setting up data statements for
c   error messages.  These names are generated automatically by PMESS,
c   the program that makes up these messages.
c LX     (*IVAA) Used for temporary storage in computing TMARKA().
c        ( formal *IVADE)  An integer array containing extra
c   information, as follows.
c  LX(1) Points to a location in Y beyond those already in use.  Values
c        of Y requested are computed at TSPECS(1) - Y(LX(1)) and stored
c        starting at Y(LX(1)+1).  If this index is 0, no more extra Y
c        values are to be computed.
c  LX(2) Index of the first equation for which the Y's above are to be
c        computed.  Y(LX(1)+1) will correspond to this first equation
c        index.
c  LX(3) Index of the last equation for which the Y's above are to be
c        computed.  Thus the Y's stored starting at Y(LX(1)+1) will
c        require no more space than half the space ordinarily required
c        for the array Y(), and may require significantly less.
c  LX(4) Maximum number of times to integrate F to get Y.  This should
c        be > 0, and less than or equal to the order of the highest
c        order differential equation.  (= 0 is allowed, but probably
c        not what you want.  It would give a value only for F.)  Space
c        must be set aside for all integrals of F, even if not all are
c        requested.  For a first order system, all Y's are just the
c        first integrals of the corresponding F's.  For higher order
c        equations, the first Y associated with a given F is the d-th
c        integral of the corresponding F, where d is the order of the
c        equation, and the last Y corresponding to the F is the first
c        integral of that F.
c  LX(5) As for LX(4), but gives the index for the fewest number of
c        times to integrate F.  Ordinarily this should be > 0.  If 0 is
c        requested, an estimate for the value of F at the delay point is
c        computed.  This should not be 0 more than once, for equations
c        covering the same index, since later such requests would write
c        over the earlier results.
c  LX(5i+k) , k = 1, 2, ... 5.  Treated as for the cases above.  If i
c        different cases of delayed Y's are to be computed, then
c        LX(5i+1) must be 0.
c LX2    (*IVADE) Value of LX(5i+2), when working on the i-th delay.
c MACT   Used in the programs which call the error message program.
c   This array difines the actions to be taken by that program.  (In
c   (*IVA,A,DA,DE,G,IN,OP)
c MACT0  (*IVADB) Used to call the message program, see MACT.
c MACT?  As for MACT, in (*IVA,CR,DB)
c MACTFV (*IVADB) As for MACT0.
c MAXDIF (*IVASC) Maximum differentiations required for stiff equations.
c MAXINT (*IVASC) Maximum integrations required.  (= max. order of
c   differential equations if equations are not stiff.)
c MAXKQ  (*IVA, BU)e
c MAXKQD (*IVAMC) Largest integration order allowed for stiff equations.
c MAXKQI (*IVAMC) Largest integ. order allowed for nonstiff equations.
c ME???? Parameters defining constants used for interaction with the
c   error message program MESS.  See comments there for definitions.
c   (In *IVA,A,DA,DE,G,IN,OP)
c METHOD (*IVAMC) Defines kind of methods being used.
c   = -1  Only stiff equations are being integrated.
c   =  0  Only nonstiff equations are being integrated.
c   =  1  Both kinds of methods are required.
c MLOC   (*IVA,A,DE) Contains locations in MTEXT for error messages.  In
c   *IVAA this data is packed using MLOC??, see below.
c MLOC?? (*IVAA) Parameters constructed to aid in making up packed data
c   for processing error messages.  Low two digits give the value of
c   KORD1I to use for the error index and later processing, the next two
c   give the error severity level, and the rest point to text used for
c   the message.
c MODF2  (*IVADB) Used in constructing the same kind of packed data as
c   described for LOCF1 above.
c MULTJ  Local to SIVAOP for calls not using F.
c MTEXT  (*IVA,A,CR,IN,OP) Text for error messages.
c MTXT?? (*IVA,A,CR,DA,DB,DE,G,IN,OP) Equivalenced into MTEXT.
c N      Used for temporary storage.  (In *IVAHC,IN,PR)
c NDTF   (*IVASC) Location in F() where difference table starts.
c NE     (*IVAMC) Number of equations in the first group.  (=NTE if
c   option 18 is not used.)
c NEDDIG (*IVADB) Parameter = -MEDDIG.
c NEPTOL (*IVAMC) Used for temporary storage and to save the value of
c   ITOLEP for error messages.
c NEQ    (formal) Total number of equations being integrated.
c NG     (*IVAMC) Used in *ivag for the number of g's in the current
c   context.
c NGSTOP (*IVAG) Dimension 2 array equivalenced to IOP6, and IOP7.  To
c   get the number of interpolating and extrapolating G-Stops.
c NGTOT  (*IVAMC) NGTOT(1) gives the number of interpolating G-Stops,
c   and NGTOT(2) gives the number of extrapolating G-Stops.
c NKDKO  (*IVASC) If this is nonzero (option 17), it gives the location
c   in KORD() where a vector defining the order of each equation is
c   specified.
c NLX    (*IVADE) Temporary index used to keep track of interpolations
c   being done to get Y() values for a delay differential equation.
c NOISEQ (*IVAMC) max(2, order of equation for which (error estimate)/
c   (error requested) is a maximum).
c NOUTKO (*IVAMC) If nonzero, gives the index in KORD where information
c   on what equations are to be included in the diagnostic output is
c   given.   See option 10.
c NSTOP  (formal) In *IVAG.  Index of the G-stop, see IFLAG.
c NTE    (*IVASC) Total number of equations being integrated = NEQ.
c NTEXT  (formao *IVADB) Character variable containing heading text.
c NTOLF  (*IVAMC) First location in F() where tolerance specifying
c   accuracy desired is stored.
c NUMDT  (*IVASC) Maximum allowed number of differences available for
c   doing an integration.
c NXTCHK (*IVA) Equivalenced to INTCHK(1), which gives the next
c   available location in INTCHK for storing data on storage allocation.
c NY     (*IVAMC) Total order of the system.
c NYNY   (*IVASC) Location in Y() where the base value for Y() is saved.
c NYNYO  (*IVADE) Equivalenced to the saved value from common of NYNY.
c OUTPUT (formal) Name of subroutine to be called for the output of
c   data or for computing G-Stops.  Not used if option 14 is set.
c OVD10  (*IVAEV) (largest floating point number) / 10.
c OVTM75 (*IVAEV) (largest floating point number) ** (-.75)
c RBQ    (*IVAMC) Array containing data for the preliminary noise test.
c RD     (formal *IVACO) Array use to contain floating point data from
c   common.  Values are returned in locations 1 to 3 as follows.
c   1    EMAX =  Max. ratio of estimated error to requested error
c   2            Reserved for future use
c   3            Reserved for future use
c REF    (*IVACR) Array of length 3 used for translating error tolerance
c   type into the factor used for exponential averaging for that type.
c RND    (*IVACR) Usually the current estimated error.  Used in deciding
c   if noise is limiting precision.
c RNOISE (*IVACR) Value used in comparison with RBQ() for preliminary
c   noise test.
c ROBND  (*IVAMC) Used to influence the selection of integration order.
c   The larger ROBND, the harder it is to increase the order and the
c   easier it is to decrease it.
c RVC2   (*IVADB) Array used for output of variables DNOISE to SNOISE in
c   common block *IVAMC.  These are variables that don't require a great
c   deal of precision.
c S      (*IVACR) Estimate of (step size) * eigenvalue of Jacobian.
c SIGMA  (*IVAMC) The k-th entry of this array contains a factor that
c   gives the amount the k-th difference is expected to increase if the
c   step size in increased.  These numbers get bigger it there is a past
c   history of increasing the step size.
c SIGMAS (*IVAA) Saved value of SIGMA(k) from the last step, where k =
c   integration order for equation with index KEMAX.
c SNOISE (*IVAMC) Value used in comparison with RBQ() on equation with
c   largest value for (error estimate) / (error requested).
c T      (formal) in *IVAIN. T(1) contains the point to be interpolated
c   to, and T(2) is used in a check that |HI| .le. |T(2)|.  When used by
c   other routines in this package, TSPECS is passed in for T.
c TB      (*IVADE) Base time for current interpolation.
c TC      (*IVADE) Original value of TN when getting past Y's for a
c   delay differential equation.
c TEMP   Used for temporary storage, in *IVAHC,PR
c TEMPA  (*IVACR) Array equivalenced to (TPS1,TPS2,TPS3,TPS4).
c TEMPAO (*IVACR) Array used to accumulate values in TEMPA.
c TG     (*IVAMC) TG(1) gives the last value of TSPECS(1) for which an
c   interpolatory G-Stop has been computed.  TG(2) is defined similarly
c   for extrapolatory G-Stops.
c TGSTOP (*IVAMC) TGSTOP(1) gives the value of TSPECS(1) where the last
c   0 for an interpolatory G-Stop was found.  TGSTOP(2) is defined
c   similarly for extrapolatory G-Stops.
c TMARK  (*IVAMC) Location of the next output point.
c TMARKA (*IVAA)  Array of length 2 equivalenced to TMARK (and TMARKX).
c TMARKX (*IVAMC) Location of the next output point to be found using
c   integration or extrapolation.  This variable must follow immediately
c   after TMARK in the common block.
c TN     (*IVASC) The value of TSPECS(1) at the conclusion of the last
c   step.
c TNEQ   (*IVADB) Array of dimension 1 equivalenced to TN so that an
c   array can be passed to *MESS.
c TOL    (formal) This is a part of F passed into *IVACR.  The first
c   location is the start of the information on the tolerances for error
c   control.
c TOLD   (*IVAG) Value of TSPECS(1) on one side of a zero.
c TOLG   (*IVAMC) Tolerance to pass to dzero when locating G-Stops.
c TOUT   (*IVAMC) Location of next output point defined by value of
c   TSPECS(3).  Such output is given with KORD(1) = 2.
c TP     (*IVA,A,DA,DE,HC) Used for temporary storage.
c TP1    (*IVAA,DA,HC,IN,PR) Used for temporary storage.
c TP2    (*IVAA,DA,HC,PR) Used for temporary storage.
c TP3    (*IVAA) Used for temporary storage.
c TPD    (*IVABU) Used for temporary storage.
c TPP    (*IVACR) Used for temporary storage.  Usually same as TPS3.
c TPS1   (*IVAA,CR) Used for temporary storage.  (In *IVACR is the
c   difference of order KQQ-2)
c TPS2   (*IVAA,CR) Used for temporary storage.  (In *IVACR is the
c   difference of order KQQ-1)
c TPS3   (*IVACR) Contains the difference of order KQQ.  This is the
c   last difference used in the corrector.
c TPS4   (*IVACR) Contains the difference of order KQQ+1.
c TPS5   (*IVACR) Temporary storage.
c TPS6   (*IVACR) Temporary storage.
c TPS7   (*IVACR) Temporary storage.
c TSAVE  (*IVAG) Value of TSPECS(1) before starting the search for a 0.
c TSPECS (formal *IVA,A,DB,DE,G)
c   TSPECS(1) is the current value of the independent variable.
c   TSPECS(2) is the current value of the step size.
c   TSPECS(3) is the increment to use between output points that give
c             output with KORD(1) = 2.
c   TSPECS(4) is the "final" output point.
c V      (*IVAMC) Array used in computing integration coefficients.
c XI     (*IVASC) XI(K) = TSPECS(1) - value of TSPECS(1) K steps
c   previous.
c W      (*IVAHC) Array used in computing integration coefficients.
c WDE    (formal, *IVADE)  Array used for working storage.  This storage
c   is used to save derivative values when iterating to get started.  To
c   be safe one should allow as much space as is allowed for differences
c   in F.  In most cases the start will not require this much space
c   however.  This array is also intended for the support of saving long
c   past histories.
c Y      (formal, *IVA,A,CR,DA,DB,DE,G,IN,PR) Array containing the
c   independent variable and all derivatives up to order one less than
c   the order of the differential equation.  Also use to save these
c   values at the beginning of the current step, the base values.
c YN     (formal, in *IVAPR)  Base values of y, these follow the
c   current values of the dependent variable, y, in Y().
c
c
c++S Default KDIM = 16
c++  Default KDIM = 20
c++  Default MAXORD = 2, MAXSTF = 1
c++  Default INTEGO, VAREQ, OUTPUT, DUMP, GSTOP, EXTRAP
c++  Default STIFF=.F., ARGM=.F., ERRSTO=.F.
c
      integer NEQ, IDIMT, IDIMY, IDIMF, IDIMK
      integer KORD(*), IOPT(*)
c--S Next line special: P=>D, X=>Q
      real             TSPECS(*), Y(*)
      real             F(*)
      external SIVAF, SIVAO
c
c *********************** Internal Variables ***************************
c
c Comments for variables used in this package can be found in the file
c   IVACOM.
c
c *********************** Type Declarations ****************************
c
      integer KDIM, MAXORD, MAXSTF
c++ Substitute for KDIM, MAXORD, MAXSTF below
      parameter (KDIM = 16, MAXORD = 2, MAXSTF = 1)
c--S Next line special: P=>D, X=>Q
      real             TN
      real             XI(KDIM)
c
c--S Next line special: P=>D, X=>Q
      real             TG(2), TGSTOP(2), TMARK, TMARKX, TOUT, TOLG
      real             ALPHA(KDIM), BETA(KDIM+1)
      real              D(MAXSTF+MAXORD,MAXORD), G(KDIM,MAXORD)
      real             V(KDIM+MAXORD)
      real             HC, HDEC, HINC, HINCC, HMAX, HMAXP9, HMIN
      real             FDAT(11)
c
      real             DS(MAXSTF+MAXORD, MAXORD), GS(KDIM)
      real             SIGMA(KDIM), RBQ(KDIM), DNOISE
      real             EAVE, EIMAX, EIMIN, EMAX, EREP, ROBND, SNOISE
c
c.    SPECIFICATION OF ENVIRONMENTAL CONSTANTS.
      real             EEPS10, EEPS16, EROV10, EEPS2
      real             EEPT75, EOVEP2, OVTM75, OVD10
      common / SIVAEV / EEPS2, EEPT75, EOVEP2, OVTM75, OVD10, EEPS10,
     1   EEPS16, EROV10
      save / SIVAEV /
      integer IOPST, KORDI, KQMAXD, KQMAXI, LDT, MAXDIF, MAXINT, NKDKO,
     1   NTE, NYNY, NDTF, NUMDT
      common / SIVASC / TN, XI, IOPST, KORDI, KQMAXD, KQMAXI, LDT,
     1   MAXDIF, MAXINT, NKDKO, NTE, NYNY, NDTF, NUMDT
c
      integer ICF,ICS,IGFLG,IGTYPE(2),IGSTOP(2),ILGREP,INGS,IOP3,IOP4,
     1   IOP5,IOP6,IOP7,IOP8,IOP9,IOP10,IOP11,IOP12,IOP13,IOP14,IOP15,
     2   IOP16,IOP17,IOP18,IOP19,IOP20,IOP21,IOP22,IOP21S,ITOLEP,IY,
     3   KEMAX,KIS,KMARK,KORD1I,KORD2I,KPRED,KQDCON,KQICON,KQMAXS,
     4   KQMXDS,KQMXIL,KQMXIP,KQMXIS,KSC,KSOUT,KSSTRT,KSTEP,LEX,LINC,
     5   LINCD,LINCQ,LSC,MAXKQD,MAXKQI,METHOD,NE,NEPTOL,NG,NGTOT,
     6   NOISEQ,NOUTKO,NTOLF,NY,IDAT(6)
      common /SIVAMC/ TG,TGSTOP,TMARK,TMARKX,TOUT,TOLG,HC,HDEC,HINC,
     1   HINCC,HMAX,HMAXP9,HMIN,ALPHA,BETA,D,G,V,DS,GS,SIGMA,RBQ,DNOISE,
     2   EAVE,EIMAX,EIMIN,EMAX,EREP,ROBND,SNOISE,FDAT,ICF,ICS,IGFLG,
     3   IGTYPE,IGSTOP,ILGREP,INGS,IOP3,IOP4,IOP5,IOP6,IOP7,IOP8,IOP9,
     4   IOP10,IOP11,IOP12,IOP13,IOP14,IOP15,IOP16,IOP17,IOP18,IOP19,
     5   IOP20,IOP21,IOP22,IOP21S,ITOLEP,IY,KEMAX,KIS,KMARK,KORD1I,
     6   KORD2I,KPRED,KQDCON,KQICON,KQMAXS,KQMXDS,KQMXIL,KQMXIP,KQMXIS,
     7   KSC,KSOUT,KSSTRT,KSTEP,LEX,LINC,LINCD,LINCQ,LSC,MAXKQD,MAXKQI,
     8   METHOD,NE,NEPTOL,NG,NGTOT,NOISEQ,NOUTKO,NTOLF,NY,IDAT
      save / SIVAMC / , / SIVASC /
c
      integer KGO, INTCHK(0:30), NXTCHK
      integer IHI, JL, J, ILOW, K, KQQ, JLIM, KFERR
c
      equivalence (INTCHK(1), NXTCHK)
      real             CM1
      parameter (CM1 = (-1.E0))
      integer IOPIVA(2)
      save IOPIVA
c
c                      Declarations for error message processing.
c
      character TEXT1(1)*10
      integer  MENTXT,MEIDAT,MEMDA1,MECONT,MERET,MEEMES,METEXT,MEIVEC
      parameter (MENTXT =23)
      parameter (MEIDAT =24)
      parameter (MEMDA1 =27)
      parameter (MECONT =50)
      parameter (MERET  =51)
      parameter (MEEMES =52)
      parameter (METEXT =53)
      parameter (MEIVEC =57)
c
      integer MACT(16), MLOC(12), MACT1(4)
c
c ********* Error message text ***************
c[Last 2 letters of Param. name]  [Text generating message.]
cAA SIVA$B
cAB The interval [1, 10**6], bounds the allowed values for NTE=$I.$E
cAC For option $I, the interval [$I, $I], bounds the allowed $C
c   values for the integration order which is set to $I.$E
cAD Option 16 must be used for error control.$E
cAE F($I) = $F, but it must be -1.0 when skipping the error check.$E
cAF For option $I, the interval [$I, $I] bounds the allowed $C
c   values for KORD($I)=$I, which is used to specify an $B
cAG output type for printing.$E
cAH output group for printing.$E
cAI equation group for variational equations.$E
cAJ order for a differential equation.$E
cAK equation group for diagnostic print.$E
cAL equation group for integration order control.$E
cAM equation group for error control.$E
cAN Option 5 argument must be .le. 0 or .gt. 4.$E
c   $
cAO KORD values for this option starting at KORD($M) are:$E
      integer LTXTAA,LTXTAB,LTXTAC,LTXTAD,LTXTAE,LTXTAF,LTXTAG,LTXTAH,
     * LTXTAI,LTXTAJ,LTXTAK,LTXTAL,LTXTAM,LTXTAN,LTXTAO
      parameter (LTXTAA=  1,LTXTAB=  7,LTXTAC= 71,LTXTAD=183,LTXTAE=226,
     * LTXTAF=290,LTXTAG=400,LTXTAH=427,LTXTAI=455,LTXTAJ=498,
     * LTXTAK=535,LTXTAL=573,LTXTAM=620,LTXTAN=655,LTXTAO=  1)
      character MTXTAA(3) * (233)
      character MTXTAB(1) * (55)
      data MTXTAA/'SIVA$BThe interval [1, 10**6], bounds the allowed val
     *ues for NTE=$I.$EFor option $I, the interval [$I, $I], bounds the$
     * allowed values for the integration order which is set to $I.$EOpt
     *ion 16 must be used for error control.$EF($I) = ','$F, but it must
     * be -1.0 when skipping the error check.$EFor option $I, the interv
     *al [$I, $I] bounds the allowed values for KORD($I)=$I, which is us
     *ed to specify an $Boutput type for printing.$Eoutput group for pri
     *nting.$Eequation gro','up for variational equations.$Eorder for a$
     * differential equation.$Eequation group for diagnostic print.$Eequ
     *ation group for integration order control.$Eequation group for err
     *or control.$EOption 5 argument must be .le. 0 or .gt. 4.$E'/
      data MTXTAB/'KORD values for this option starting at KORD($M) are:
     *$E'/

c End of automatically generated error message code.
c
c        for KGO =     1      2      3      4      5      6      7
      data MLOC / LTXTAI,LTXTAK,LTXTAL,LTXTAM,LTXTAJ,LTXTAG,LTXTAH,
     2            LTXTAB,LTXTAC,LTXTAD,LTXTAE, LTXTAN /
c           KGO        8      9     10     11      12
c
c                      1  2  3 4       5  6       7       8       9 10
      data MACT / MEEMES,38,24,0, MENTXT, 0, METEXT, MECONT, MEMDA1,0,
     1  METEXT, MEIDAT,0, MEIVEC,0, MECONT /
c           11      12 13     14 15     16
      data MACT1 / METEXT, MEIVEC, 0, MERET /
      data IOPIVA(1) / 1111 /
      data TEXT1 / 'IOPT()= $B' /
c
c ************** START OF EXECUTABLE CODE ******************
c
c     **** TEST IF CONTINUING AN INTEGRATION
      if (KORD(1) .ne. 0) go to 330
c     **** INITIALIZE VARIOUS SCALARS
      KSTEP = 0
      KQMXIS = 0
      KORD2I = -5
      KORD(2) = -1
      NTE = NEQ
      NE = NTE
      TOLG = 0.E0
c     **** SET UP OPTIONS
      if (IOPT(1) .ne. 0) call SIVAOP(IOPT, F)
      call SIVAOP(IOPIVA, F)
      if (IOPT(1) .eq. 0) IOPIVA(2) = 1
c
      if ((NE .le. 0) .or. (NE .gt. 1000000)) then
         IDAT(1) = NE
         KGO = 8
         go to 650
      end if
c                         Set up diagnostic print on storage allocation.
      INTCHK(0) = 245
      if (IOP10 .ne. 0) INTCHK(0) = 247
c
c     **** CHECK TSPECS STORAGE ALLOCATION
      INTCHK(2) = IDIMT
      INTCHK(3) = 4
      NXTCHK = 4
      if (IOP5 .ne. 0) then
         INTCHK(4) = 5
         INTCHK(5) = 5
         if (IOP5 .gt. 0) then
            INTCHK(6) = IOP5 - 4
            if (IOP5 .lt. 5) then
               KGO = 12
               go to 600
            end if
         else
            IHI = -IOP5
            JL = 4
            do 15 IHI = IHI, IDIMK-3, 3
               J = abs(KORD(IHI))
               if (J .eq. 0) go to 20
               if (abs(KORD(IHI + 2)) .gt. 1) then
                  IDAT(2) = -1
                  IDAT(3) = 1
                  KGO = 6
                  go to 600
               end if
               if ((J .le. JL) .or. (J .gt. KORD(IHI+1))) then
                  KGO = 7
                  IDAT(2) = JL + 1
                  IDAT(3) = KORD(IHI+1)
                  go to 610
               end if
               JL = KORD(IHI+1)
   15       continue
            if (KORD(IHI) .ne. 0) IHI = IHI + 3
   20       INTCHK(6) = JL - 4
         end if
         NXTCHK = 7
      end if
   25 call OPTCHK(INTCHK, IOPT, 'SIVA / TSPECS$E')
      if (NXTCHK .lt. 0) KORD2I = -4
c
c     **** CHECK KORD STORAGE ALLOCATION
      INTCHK(2) = IDIMK
      INTCHK(3) = NE + 3
      NXTCHK = 4
      if (IOP5 .lt. 0) then
         INTCHK(4) = 5
         INTCHK(5) = -IOP5
         INTCHK(6) = IHI + IOP5
         NXTCHK = 7
      end if
c
c++  Code for VAREQ is active
      if (IOP18 .ne. 0) then
         NE = abs(KORD(IOP18))
         INTCHK(NXTCHK) = 18
         ILOW = IOP18
         KGO = 1
c.       **** CHECK OPTION FOR VALID INPUT
         go to 430
      end if
c++  End
   30 continue
      if (NKDKO .ne. 0) then
c                        **** STORAGE ALLOCATED FOR ODE ORDERS
         INTCHK(NXTCHK) = 17
         INTCHK(NXTCHK+1) = NKDKO
         INTCHK(NXTCHK+2) = NTE
         NXTCHK = NXTCHK + 3
      end if
c++  Code for STIFF is inactive
c      IF (IOPST .ne. 0) then
c         INTCHK(NXTCHK) = 17
c         INTCHK(NXTCHK+1) = IOPST
c         INTCHK(NXTCHK+2) = NTE
c         NXTCHK = NXTCHK + 3
c      end if
c++  End
c
c **** SET INITIAL INTEGRATION ORDERS, TEST ODE ORDERS ****
c
      MAXINT = 0
      MAXDIF = 0
      NY = 0
      do 80 K = 1, NTE
         if (NKDKO .ne. 0) KORDI = KORD(NKDKO + K - 1)
         NY = NY + abs(KORDI)
c++  Code for STIFF is inactive
c      IF (IOPST .EQ. 0) GO TO 60
cc.    **** CHECK FOR POTENTIAL STIFF EQUATION
c      JS = abs(KORD(IOPST+K-1)) - 1
c      IF ( JS ) 52,60,54
cc.    **** EQUATION IS NOT ACTIVE
c   52 KQQ = 0
c      GO TO 56
cc.    **** EQUATION USES IMPLICIT METHOD
c   54 KQQ = -1
c      IF (JS .GT. abs(KORDI)) then
c        Set up an error message.
c      end if
c      MAXINT = max(MAXINT, abs(KORDI) - JS)
c   56 IF (KORDI .GE. 0) GO TO 70
c      KORDI = -1 - KORDI
c      JS = JS - 1
c      MAXDIF = max(MAXDIF, JS, 1)
c      GO TO 70
c++  End
c     **** EQUATION IS TO USE AN EXPLICIT METHOD
   60    KQQ = 1
         MAXINT = max(MAXINT, KORDI)
   70    if ((KORDI .gt. MAXORD) .or. (KORDI .le. 0)) then
c                    Set up error message.  KORDI is out of range.
            IDAT(1) = 17
            IDAT(2) = 1
            IDAT(3) = MAXORD
            if (NKDKO .ne. 0) then
               KGO = 5
               ILOW = NKDKO
               IHI = NKDKO + K - 1
               go to 640
            else
               KGO = 9
               IDAT(4) = KORDI
               go to 650
            end if
         end if
         KORD(K + 3) = KQQ
   80 continue
c     **** SET FLAGS WHICH DEPEND ON METHOD USED
c++  Code for STIFF is inactive
c      METHOD = 1
c      IF (MAXINT .GT. 0) IF (MAXDIF) 85,90,85
c      METHOD = -1
c   85 CONTINUE
c      KPRED = 5
c      GO TO 100
c++  End
   90 METHOD = 0
      KPRED = 1
  100 continue
c
c ******* CHECK KORD FOR DIAGNOSTIC OUTPUT CONTROL *********
c
c++  Code for OUTPUT is active
      if (IOP10 .gt. 0) then
         if (NOUTKO .ne. 0) then
            INTCHK(NXTCHK) = 10
            ILOW = NOUTKO
c.    **** Check option for valid input
            KGO = 2
            go to 430
         end if
      end if
c++  End
  110 continue
c
c ********** CHECK KORD FOR INTEGRATION ORDER CONTROL ******
c
c++  Code for INTEGO is active
      if (IOP19 .ne. 0) then
c.           **** Check option for valid input
         INTCHK(NXTCHK) = 19
         ILOW = IOP19
         JLIM = -30
         KGO = 3
         go to 430
      end if
c++  End
  120 continue
c
c ********** CHECK SET UP FOR ERROR TOLERANCES *************
c
      INTCHK(NXTCHK) = 16
      ILOW = IOP16
      JLIM = -5
      KGO = 4
      if (IOP16 .ne. 0) go to 430
c.                      **** IN CURRENT CODE, IOP16=0 IS AN ERROR
      KGO = 10
      go to 650
  150 continue
c     **** CHECK KORD STORAGE ALLOCATION
      call OPTCHK(INTCHK, IOPT, 'SIVA / KORD$E')
      if (NXTCHK .lt. 0) KORD2I = -4
c
c     ******** DONE CHECKING KORD STORAGE ALLOCATION *******
c
c     **** CHECK  Y  STORAGE ALLOCATION
      INTCHK(2) = IDIMY
      INTCHK(3) = NY + NY
      NXTCHK = 4
      NYNY = NY + 1
      call OPTCHK(INTCHK, IOPT, 'SIVA / Y$E')
      if (NXTCHK .lt. 0) KORD2I = -4
c
c     **** CHECK  F  STORAGE ALLOCATION
      INTCHK(2) = IDIMF
      INTCHK(3) = NTE
      NXTCHK = 4
      if (IOP16 .ne. 0) then
c                                Error tolerance info.
         INTCHK(4) = 16
         INTCHK(5) = NTOLF
         INTCHK(6) = IHI - IOP16 + 1
         NXTCHK = 7
      end if
      if (IOP12 .gt. 0) then
         INTCHK(NXTCHK) = 12
         INTCHK(NXTCHK+1) = IOP12
         INTCHK(NXTCHK+2) = 4
         NXTCHK = NXTCHK + 3
      end if
      if (IOP21 .gt. 0) then
         INTCHK(NXTCHK) = 21
         INTCHK(NXTCHK+1) = IOP21
         INTCHK(NXTCHK+2) = 1
         NXTCHK = NXTCHK + 3
      end if
c
c++  Code for ERRSTO is inactive
c      IF (IOP20 .ne. 0) then
cc.                                Space for saving error estimates
c         INTCHK(NXTCHK) = 20
c         INTCHK(NXTCHK) = IOP20
c         INTCHK(NXTCHK) = NTE
c         NXTCHK = NXTCHK + 3
c      end if
c++  Code for STIFF is inactive
c      if (IOP21 .gt. 0) then
cc.                               Info. for stiff equations
c         INTCHK(NXTCHK) = 21
c         INTCHK(NXTCHK+1) = IOP21
c         INTCHK(NXTCHK+2) = IOP21S
c         NXTCHK = NXTCHK + 3
c      end if
c      MAXKQD = min(MAXKQI, 6)
c++  End
c                          Set aside space for the difference tables.
      INTCHK(NXTCHK) = 0
      INTCHK(NXTCHK+1) = -KDIM * NTE
      INTCHK(NXTCHK+2) = 0
      NXTCHK = NXTCHK + 3
      INTCHK(NXTCHK) = -5 * NTE
      call OPTCHK(INTCHK, IOPT, 'SIVA / F$E')
      if (NXTCHK .lt. 0) then
         KORD2I = -4
      else if (KORD2I .ne. -4) then
         do 290 K = NXTCHK+1, INTCHK(NXTCHK)
            if (INTCHK(INTCHK(K)) .eq. 0) then
               NDTF = INTCHK(INTCHK(K)+1)
               NUMDT = min(KDIM, (INTCHK(INTCHK(K)+2)-NDTF+1) / NTE)
               MAXKQI = NUMDT - 1
            else
c         Take a quick return if needed space was not specified by user.
               KORD2I = -4
            end if
  290    continue
      end if
      if (IOP9 + abs(IOP10) + IOP11 .ne. 0) then
c Insure user doesn't get in trouble with F not iniitalized.
         do 300 K = NDTF, NDTF + NTE*NUMDT - 1
            F(K) = 0.E0
  300    continue
      end if
  320 continue
      if ((KORD2I .eq. -4) .or. (IOP10 .ne. 0)) then
         MACT1(3) = IOPIVA(2)
         call MESS(MACT1, TEXT1, IOPT)
         KORD1I = 24
         KORD(1) = 24
      end if
      TMARK = TSPECS(1)
      TMARKX = TSPECS(4) + TSPECS(2)
c
c     **** DONE WITH INITIALIZATION AND CHECKING INPUTS
      if (IOP13 + IOP14 + IOP15 .ne. 0) return
  330 call SIVAA(TSPECS, Y, F, KORD, SIVAF, SIVAO)
      return
c
c ************ LOOP TO CHECK OPTION SPECIFICATIONS *********
c
  430 JL = 0
      do 560 IHI = ILOW, IDIMK
         J = KORD(IHI)
         go to (460, 480, 490, 490), KGO
c     **** CHECK ON VARIATIONAL EQUATIONS
  460    continue
c++  Code for VAREQ is active
         if (J - NTE) 470, 565, 620
  470    if (J .eq. 0) go to 560
         if (J .le. JL) go to 620
c++  End
c     **** Check on diagnostic output option
  480    continue
c++  Code for OUTPUT is active
c.    **** CHECK IF DONE
         if (J .ge. NTE) go to 565
         if (J .le. JL) go to 620
         go to 550
c++  End
  490    continue
c     **** Check integration order control (KGO=3) and
c     **** error tolerance equation grouping (KGO=4).
         if (J - NTE) 500, 565, 620
  500    if (J) 510, 530, 540
  510    if ((JL .le. 0) .and. (IHI .ne. ILOW)) go to 620
         if (J .lt. JLIM) then
c                         Output an error message.
            IDAT(2) = JLIM
            IDAT(3) = 0
            go to 630
         end if
  520    JL = -JL
         go to 560
  530    if (KGO .eq. 3) go to 520
         KFERR = NTOLF + IHI - ILOW
         if (F(KFERR) .eq. CM1) go to 510
c                         Set up error message, TOL must be -1.
            IDAT(1) = KFERR
            KGO = 11
            go to 650
  540    if (abs(JL) .ge. abs(J)) go to 620
  550    JL = J
  560    continue
  565 NXTCHK = NXTCHK + 3
      INTCHK(NXTCHK-2) = ILOW
      INTCHK(NXTCHK-1) = IHI - ILOW + 1
      go to (30, 110, 120, 150), KGO
c
c     **** AN ERROR HAS BEEN MADE
c                  Error in setting up TSPECS for extra output
  600 IHI = IHI + 2
  610 ILOW = -IOP5
      go to 630
c                  Error in KORD indices
  620 IDAT(2) = abs(JL) + 1
      IDAT(3) = NTE
c                  Set up for print of message about KORD
  630 IDAT(1) = INTCHK(NXTCHK)
  640 IDAT(4) = IHI
      IDAT(5) = KORD(IHI)
c
c ***************** Process Errors *************************************
c
  650 KORD2I = -4
      MACT(4) = LTXTAF
      if (KGO .ge. 8) MACT(4) = -1
      MACT(6) = MLOC(KGO)
c--S Next line special: P=>S, X=>D
      CALL SMESS(MACT, MTXTAA, IDAT, FDAT)
      if (KGO .lt. 8) then
         MACT(10) = ILOW
         MACT(13) = ILOW
         MACT(15) = -min(IHI+2, IDIMK)
         CALL MESS(MACT(9), MTXTAB, KORD)
         if (KGO .le. 4) go to 565
      end if
c              5   6   7    8    9   10   11  12
      go to (100, 25, 25, 320, 100, 150, 660, 25), KGO - 4
  660 KGO = 4
      go to 565
      end
c   End of SIVA

      subroutine SIVAA(TSPECS, Y, F, KORD, SIVAF, SIVAO)
c>> 1989-02-24 SIVAA  Krogh   Big error with BETA(2)=1+epsilon -- looped
c>> 1988-07-27 SIVAA  Krogh   Fixed to allow restart on a restart.
c>> 1988-03-07 SIVAA  Krogh   Initial code.
c
c  MAIN SUBROUTINE FOR VARIABLE ORDER INTEGRATION OF ORDINARY
c  DIFFERENTIAL EQUATIONS
c
      integer KORD(*)
c--S Next line special: P=>D, X=>Q
      real             TSPECS(*), Y(*)
      real             F(*)
      external SIVAF, SIVAO
c
      integer KDIM, MAXORD, MAXSTF
c++ Substitute for KDIM, MAXORD, MAXSTF below
      parameter (KDIM = 16, MAXORD = 2, MAXSTF = 1)
c--S Next line special: P=>D, X=>Q
      real             TN
      real             XI(KDIM)
c
c--S Next line special: P=>D, X=>Q
      real             TG(2), TGSTOP(2), TMARK, TMARKX, TOUT, TOLG
      real             ALPHA(KDIM), BETA(KDIM+1)
      real              D(MAXSTF+MAXORD,MAXORD), G(KDIM,MAXORD)
      real             V(KDIM+MAXORD)
      real             HC, HDEC, HINC, HINCC, HMAX, HMAXP9, HMIN
      real             FDAT(11)
c
      real             DS(MAXSTF+MAXORD, MAXORD), GS(KDIM)
      real             SIGMA(KDIM), RBQ(KDIM), DNOISE
      real             EAVE, EIMAX, EIMIN, EMAX, EREP, ROBND, SNOISE
c
c.    SPECIFICATION OF ENVIRONMENTAL CONSTANTS.
      real             EEPS10, EEPS16, EROV10, EEPS2
      real             EEPT75, EOVEP2, OVTM75, OVD10
      common / SIVAEV / EEPS2, EEPT75, EOVEP2, OVTM75, OVD10, EEPS10,
     1   EEPS16, EROV10
      save / SIVAEV /
      integer IOPST, KORDI, KQMAXD, KQMAXI, LDT, MAXDIF, MAXINT, NKDKO,
     1   NTE, NYNY, NDTF, NUMDT
      common / SIVASC / TN, XI, IOPST, KORDI, KQMAXD, KQMAXI, LDT,
     1   MAXDIF, MAXINT, NKDKO, NTE, NYNY, NDTF, NUMDT
c
      integer ICF,ICS,IGFLG,IGTYPE(2),IGSTOP(2),ILGREP,INGS,IOP3,IOP4,
     1   IOP5,IOP6,IOP7,IOP8,IOP9,IOP10,IOP11,IOP12,IOP13,IOP14,IOP15,
     2   IOP16,IOP17,IOP18,IOP19,IOP20,IOP21,IOP22,IOP21S,ITOLEP,IY,
     3   KEMAX,KIS,KMARK,KORD1I,KORD2I,KPRED,KQDCON,KQICON,KQMAXS,
     4   KQMXDS,KQMXIL,KQMXIP,KQMXIS,KSC,KSOUT,KSSTRT,KSTEP,LEX,LINC,
     5   LINCD,LINCQ,LSC,MAXKQD,MAXKQI,METHOD,NE,NEPTOL,NG,NGTOT,
     6   NOISEQ,NOUTKO,NTOLF,NY,IDAT(6)
      common /SIVAMC/ TG,TGSTOP,TMARK,TMARKX,TOUT,TOLG,HC,HDEC,HINC,
     1   HINCC,HMAX,HMAXP9,HMIN,ALPHA,BETA,D,G,V,DS,GS,SIGMA,RBQ,DNOISE,
     2   EAVE,EIMAX,EIMIN,EMAX,EREP,ROBND,SNOISE,FDAT,ICF,ICS,IGFLG,
     3   IGTYPE,IGSTOP,ILGREP,INGS,IOP3,IOP4,IOP5,IOP6,IOP7,IOP8,IOP9,
     4   IOP10,IOP11,IOP12,IOP13,IOP14,IOP15,IOP16,IOP17,IOP18,IOP19,
     5   IOP20,IOP21,IOP22,IOP21S,ITOLEP,IY,KEMAX,KIS,KMARK,KORD1I,
     6   KORD2I,KPRED,KQDCON,KQICON,KQMAXS,KQMXDS,KQMXIL,KQMXIP,KQMXIS,
     7   KSC,KSOUT,KSSTRT,KSTEP,LEX,LINC,LINCD,LINCQ,LSC,MAXKQD,MAXKQI,
     8   METHOD,NE,NEPTOL,NG,NGTOT,NOISEQ,NOUTKO,NTOLF,NY,IDAT
      save / SIVAMC / , / SIVASC /
c
      real             CMP75, C0, C1M5, C1M3, C2P5M3, C8M3, CP0625, CP1
      real             CP25, CP3, CP4, CP5, CP875, C1, C1P125, C1P3, C2
      real             C6, C10, C16, C4096
      parameter (CMP75 = -.75E0)
      parameter (C0 = 0.E0)
      parameter (C1M5 = 1.E-5)
      parameter (C1M3 = 1.E-3)
      parameter (C2P5M3 = 2.5E-3)
      parameter (C8M3 = 8.E-3)
      parameter (CP0625 = .0625E0)
      parameter (CP1 = .1E0)
      parameter (CP25 = .25E0)
      parameter (CP3 = .3E0)
      parameter (CP4 = .4E0)
      parameter (CP5 = .5E0)
      parameter (CP875 = .875E0)
      parameter (C1 = 1.E0)
      parameter (C1P125 = 1.125E0)
      parameter (C1P3 = 1.3E0)
      parameter (C2 = 2.E0)
      parameter (C6 = 6.E0)
      parameter (C10 = 10.E0)
      parameter (C16 = 16.E0)
      parameter (C4096 = 4096.E0)
c
      integer LDIS, KEXIT, I, J, K, J1, J2, L, LX
      real             TP, TP1, TP2, TP3, HH, DISADJ
      real             SIGMAS, TPS1, TPS2, EIMINO, EXR
c--S Next line special: P=>D, X=>Q
      real              TMARKA(2), XP, XP1
      equivalence (G(1, 1), HH), (TMARKA(1), TMARK)
      equivalence (KEXIT, IOP17)
      save EIMINO, EXR, SIGMAS, TP, TP1, TP2, TP3, TPS1, TPS2, DISADJ,
     1   LDIS
c
c                      Declarations for error message processing.
c
      integer MACT(17), MLOC(8), MENTXT, MERET, MEEMES, METEXT
      parameter (MENTXT =23)
      parameter (MERET  =51)
      parameter (MEEMES =52)
      parameter (METEXT =53)
c
c ********* Error message text ***************
c[Last 2 letters of Param. name]  [Text generating message.]
cAA SIVAA$B
cAB At: TN=$F, KSTEP=$I, with H=$F$E
cAC A previously reported error was fatal.$E
cAD Print points not properly ordered: TSPEC($I)=$F$E
cAE An error tolerance of 0 requires setting special flags.  $B
cAF Step size reduced too fast, doing a restart.  $B
cAG H is so small that TN + H = TN.  $B
cAH Error tolerance too small.  $B
cAI Step size at end of start < HMIN=$F, $B
cAJ Error estimates require a stepsize < HMIN=$F, $B
cAK (Estimated Error) / (Requested Error) for equation $I is $F.  $B
cAL Tolerance $I is F($I) = $F.$B
cAM Tolerance $I is F($I) * F($I) = $F * $F = $F.$B
cAN   Replacing F($I) with $F.$E
      integer LTXTAA,LTXTAB,LTXTAC,LTXTAD,LTXTAE,LTXTAF,LTXTAG,LTXTAH,
     * LTXTAI,LTXTAJ,LTXTAK,LTXTAL,LTXTAM,LTXTAN
      parameter (LTXTAA=  1,LTXTAB=  8,LTXTAC= 40,LTXTAD= 80,LTXTAE=129,
     * LTXTAF=189,LTXTAG=237,LTXTAH=272,LTXTAI=302,LTXTAJ=342,
     * LTXTAK=390,LTXTAL=455,LTXTAM=484,LTXTAN=531)
      character MTXTAA(3) * (186)
c
      integer LOCM, MLOCAC, MLOCAD, MLOCAE, MLOCAF, MLOCAG, MLOCAH,
     1   MLOCAI, MLOCAJ

      parameter (LOCM = 32 * 256)
c                     KORD1I   Severity   Loc. message
      parameter (MLOCAC = 23 + 32 * (99 + 256 * LTXTAC))
      parameter (MLOCAD = 13 + 32 * (38 + 256 * LTXTAD))
      parameter (MLOCAE = 22 + 32 * (38 + 256 * LTXTAE))
      parameter (MLOCAF =  3 + 32 * (14 + 256 * LTXTAF))
      parameter (MLOCAG = 21 + 32 * (38 + 256 * LTXTAG))
      parameter (MLOCAH =  2 + 32 * (25 + 256 * LTXTAH))
      parameter (MLOCAI = 11 + 32 * (38 + 256 * LTXTAI))
      parameter (MLOCAJ = 12 + 32 * (38 + 256 * LTXTAJ))
c
      data MTXTAA/'SIVAA$BAt: TN=$F, KSTEP=$I, with H=$F$EA previously r
     *eported error was fatal.$EPrint points not properly ordered: TSPEC
     *($I)=$F$EAn error tolerance of 0 requires setting special flags. $
     * ','$BStep size reduced too fast, doing a restart.  $BH is so smal
     *l that TN + H = TN.  $BError tolerance too small.  $BStep size at$
     * end of start < HMIN=$F, $BError estimates require a steps','ize <
     * HMIN=$F, $B(Estimated Error) / (Requested Error) for equation $I$
     * is $F.  $BTolerance $I is F($I) = $F.$BTolerance $I is F($I) * F(
     *$I) = $F * $F = $F.$B  Replacing F($I) with $F.$E'/
      data MLOC / MLOCAC, MLOCAD, MLOCAE, MLOCAF, MLOCAG, MLOCAH,
     1   MLOCAI, MLOCAJ /
c
c                      1 2 3 4       5 6       7       8 9      10
      data MACT / MEEMES,0,0,0, MENTXT,0, METEXT, MENTXT,0, METEXT,
     1  MENTXT, 0, METEXT, MENTXT, LTXTAN, METEXT, MERET /
c           11 12      13      14      15      16     17
c
      data EXR, EIMINO / CP1, C8M3 /
      data LDIS / 0 /
c
c ************** START OF EXECUTABLE CODE ******************************
c
  660 if (KORD2I) 670, 1380, 1840
  670 if (KORD2I .eq. -1) go to 2140
      if (KORD2I .eq. -5) go to 720
      if (KORD2I+8) 1840, 710, 1840
c     **** SPECIAL OUTPUT CASE (EXTRAPOLATION OR GSTOP)
  680 go to (1220, 1190, 830, 690, 690, 2190), KEXIT
c     **** RESET TMARK BEFORE GOING WHERE DIRECTED BY KEXIT
  690 KEXIT = KEXIT - 2
      go to 1890
c     **** USER HAS REQUESTED A RESTART
  700 KORD2I = -5
      IGFLG = 0
      if ((KORD(2) .ge. 2) .and. (LSC .lt. 3)) KORD2I = -8
      if ((KORD1I .le. 3) .or. (KORD1I .eq. 5)) go to 1890
      if (KORD2I .eq. -5) go to 720
c                Set up for a discontinuity
  710 XP = TSPECS(1)
      if (KORD(2) .eq. 3) then
c                Adjust for discontinuity in Y
         J = 0
         do 712 I = 1, NTE
            if (NKDKO .ne. 0) KORDI = KORD(NKDKO + I -1)
            K = 1
            J = J + KORDI
            XP1 = Y(J)
  711       Y(NYNY + J - K) = Y(NYNY + J - K) + XP1
            if (K .lt. KORDI) then
               XP1 = Y(J-K) + (TN - XP) * XP1 / real(K)
               K = K + 1
               go to 711
            end if
  712    continue
      else
         IGFLG = -3
      end if
      DISADJ = HH
      XP1 = (XP - TN) / XI(1)
      if (XP1 .lt. CP25) then
         K = 1
         if (XP1 .lt. CMP75) K = 2
         TSPECS(1) = TN - XI(K)
         TSPECS(2) = 2.E0*(TN - TSPECS(1))
         call SIVAIN(TSPECS(1), Y, F, KORD)
         do 713 J = 1, NY
            Y(NYNY + J - 1) = Y(J)
  713    continue
c          Move difference tables back one step
         TN = TSPECS(1)
  714    call SIVABU(F, KORD)
         if (K .eq. 2) then
            KSC = max(KSC - 1, 1)
            do 715 K = max(1, KSC), IOP11-1
               BETA(K+1) = BETA(K) * (XI(K) / (XI(K+1) - XI(1)))
  715       continue
            K = 0
            go to 714
         end if
      end if
c      Take step to discontinuity.
      HH = (XP - TN)
      EREP = -abs(EREP)
      KIS = 1000
      LDIS = 1
      LSC = 0
      LINC = 0
      HINCC = C1P125
      KSSTRT = KSTEP + 2
      go to 1115
c ********
c INITIALIZE FOR STARTING AN INTEGRATION
c ********
  720 HH = TSPECS(2)
      LSC = 8
      LINC = -5
      LINCD = 64
      LDT = -4
      KQMAXI = 0
      EIMIN = CP3
      EAVE = C10
      XI(1) = C0
      IGFLG = 0
      KIS = 0
      KSSTRT = 0
      ROBND = 0.E0

c GO COMPUTE INITIAL DERIVATIVES
  730 KORD2I = -6
c++  Code for VAREQ is active
      ICS = 0
      ICF = NE
c++  End
  740 KORD1I = KPRED
      go to 1360
c   RETURN AFTER COMPUTING INITIAL (OR NOISE TEST) DERIVATIVES
  750 continue
c++  Code for VAREQ is active
      if (IOP18 .eq. 0) go to 790
c.    **** SPECIAL LOGIC TO COMPUTE VARIATIONAL DERIVATIVES
  760 if (KORD1I .eq. 3) go to 770
      KORD1I = 3
      KORD(3) = 0
  770 if (ICF .eq. NTE) go to 790
      if (KORD2I .eq. -6) then
         if (ICS .eq. 1) go to 790
      end if
      ICS = ICF + 1
  780 KORD(3) = KORD(3) + 1
      ICF = IOP18 + KORD(3)
      ICF = abs(KORD(ICF))
      if (ICF) 1360, 780, 1360
c++  End
  790 ICS = 1
      ICF = NE
c++  Code for VAREQ is active
      if (KORD2I .eq. 0) if (EREP) 2220, 2220, 1430
c++  End
      if (LINC + 5) 1490, 810, 1490
c END OF SPECIAL CODE FOR INITIALIZATION AT THE START
c ********
c UPDATE VARIABLES TO PREPARE FOR NEXT STEP
c ********
  800 LDT = 0
      EIMIN = C2P5M3+EIMIN*(C6*EAVE+EIMAX)/((EIMIN+C1)*(C6*EIMAX+EAVE))
  810 do 820 J = 1, NY
  820    Y(NYNY + J - 1) = Y(J)
      TN = TSPECS(1)
c ********
c TEST FOR VARIOUS TYPES OF OUTPUT
c ********
c     TEST FOR END OF STEP OUTPUT (OR IF DUMP OUTPUT TO BE TESTED FOR)
      if (IOP6 .eq. 0) go to 840
c   SET UP FOR END OF STEP OUTPUT
      KORD1I = 6
      go to 1810
c     SET UP AFTER OTHER COMPUTATIONS HAVE BEEN MADE (TSPECS(1).NE.TN)
  830 KORD1I = 6
      go to 860
c     TEST FOR AN OUTPUT POINT
  840 if (HH * (TMARK - TN)) 1760, 1760, 850
c     TEST FOR TOO MANY STEPS OUTPUT
  850 continue
      if (KSOUT .gt. KSTEP) go to 890
      KORD1I = 4
  860 if (TSPECS(1) .eq. TN) go to 1810
c     GO INTERPOLATE VALUES AT END OF LAST STEP
  870 TSPECS(1) = TN
      go to 1780
c     CONTINUE AFTER TOO MANY STEPS OUTPUT
  880 continue
      KSOUT = KSTEP + IOP4
  890 continue
c++  Code for DUMP is active
      if (IOP9 .eq. 0) go to 920
c++  Code for DUMP & STIFF is inactive
c      KQMXDS=KQMAXD
c++  Code for DUMP is active
      KQMXIS = KQMAXI
      KIS = KIS + 1
      if (KIS .ne. 0) go to 920
c.   TIME TO DUMP THE SOLUTION
  900 KORD1I = 9
c.    SET TO UPDATE THE DIFFERENCE TABLE
      if (LDT .eq. 1) go to 1810
c Note that doing this update can lead to very small differences in the
c results because of round off differences.
      LDT = -2
      go to 1780
c++  End
c RETURN AFTER DUMPING THE SOLUTION
  910 continue
c++  Code for DUMP is active
      KIS = 2
c.    TEST IF SOLUTION DUMP DUE TO RESTART, END, OR
c.    DROP IN INTEG. ORDER
      if (LINC .lt. 0) if (LINC + 6) 1860, 1750, 1180
c++  End
c END OF TESTING FOR VARIOUS TYPES OF OUTPUT
c   TEST IF STEPSIZE MAY BE INCREASED OR IF TESTS SHOULD BE MADE
c   FOR DECREASING THE STEPSIZE (OR IF STARTING OR USER SELECTING H)
  920 KSTEP = KSTEP + 1
      if (LINC) 930, 980, 940
  930 continue
c     **** ONLY POSSIBLE VALUES AT THIS POINT ARE
c          LINC = -2 OR -5
      if (LINC + 5) 1120, 1110, 1120
c ********
c ERROR ESTIMATES INDICATE STEPSIZE CAN BE INCREASED
c ********
  940 HC = HINCC
      if (LINC .gt. 1) HC = HC ** LINC
      HH = HC * HH
c     TEST IF NEW STEPSIZE IS TOO BIG
      if (abs(HH) .gt. HMAX) if (HMAX) 970, 970, 960
  950 EAVE = EIMAX
      ROBND = CP3 + HINCC
      go to 1110
c     NEW STEPSIZE IS TOO BIG
  960 if (abs(XI(1)) .ge. HMAXP9) go to 970
      HH = sign(HMAX, HH)
      go to 950
c     RESTORE THE OLD STEPSIZE
  970 HH = XI(1)
      LINC = 0
      go to 1150
c END OF CODE FOR CASE WHEN ERROR ESTIMATES INDICATE STEPSIZE INCREASE
c ********
c TEST IF ESTIMATED ERRORS INDICATE STEPSIZE SHOULD BE DECREASED
c ********
  980 ROBND = C1P3
      if (EIMAX .le. EAVE) go to 990
      EAVE = EAVE + CP4 * (EIMAX - EAVE)
      if ((EIMAX * EMAX) - C1M3) 1000, 1010, 1010
  990 EAVE = EIMAX
      if ((EIMAX * EMAX) .ge. EIMIN) go to 1010
 1000 ROBND = CP3 + (SIGMA(KQMAXS)/SIGMA(KQMAXS-1))
      go to 1180
c     TEST IF STEPSIZE SHOULD BE REDUCED
 1010 if (EMAX * EIMAX .lt. EXR * EAVE) go to 1180
c ********
c ERROR ESTIMATES INDICATE STEPSIZE SHOULD BE REDUCED
c ********
      HC = HDEC
      if (EIMIN .le. EIMINO) go to 1030
      EIMIN = EIMINO
      HC = CP875
 1030 HH = HC * XI(1)
      if (LSC - 1) 1040, 1080, 1090
 1040 if (abs(HH) .ge. HMIN) go to 1090
      if (abs(CP875*XI(1)) .le. HMIN) if (LINC) 1050, 970, 970
      HH = sign(HMIN, HH)
      go to 1090
c     STEPSIZE IS TOO SMALL TO BE REDUCED
c     SET UP ERROR INDICATORS AND PREPARE FOR RETURN TO USER
 1050 KORD1I = 8
      go to 2240
c     PROCEED WITH CURRENT STEPSIZE DESPITE ERROR BEING TOO BIG
 1070 HH = XI(1)
      TSPECS(2) = HH
      EMAX = C0
      LINC = 0
      if (KORD1I - 2) 1420, 1430, 1430
c     SET LSC TO END STARTING PHASE
 1080 LSC = 2
c     CHECK IF REPEATING A STEP
 1090 if (LINC .ne. -1) go to 1110
c   WHEN REPEATING A STEP, BACK UP THE DIFFERENCES AND STEPSIZE INFO.
 1100 call SIVABU(F, KORD)
c   TEST IF NOISE TEST (LINC = -7) OR IF H IS NOT
c     BEING CHANGED (LINC = -4)
      if (LINC + 4) 1780, 1180, 1110
c ********
c STEPSIZE IS TO BE CHANGED
c ********
 1110 continue
c MODIFY STEPSIZE TO REDUCE ROUNDOFF ERROR IN ACCUMULATING INDEP. VAR.
      TP = C2 * abs(TN) + C4096 * abs(HH)
      TP = (TP + abs(HH)) - TP
      if (TP .ne. C0) HH = sign(TP, HH)
c     TEST IF NEW STEPSIZE SELECTED ACTUALLY GIVES A CHANGE
      if (HH .eq. TSPECS(2)) go to 1140
 1115 TSPECS(2) = HH
      if (IOP8 .eq. 0) go to 1140
c     SETUP TO TELL USER ABOUT STEPSIZE CHANGE (OR TO CHANGE STEPSIZE)
 1120 KORD1I = 8
      go to 860
c     RETURN AFTER TELLING USER ABOUT STEPSIZE CHANGE
 1130 HH = TSPECS(2)
 1140 if (HH .ne. XI(1)) KQICON = -1
 1150 HC = min(EOVEP2, abs(HH)) / EEPS2
c ********
c PREPARE FOR BEGINNING A NEW STEP
c ********
      if (LINC .gt. 0) then
         LINC = min(LINC, LINCQ) + LINCQ
         go to 1190
      end if
 1180 LINC = LINCD
 1190 if (HC .gt. abs(TN)) go to 1200
c     **** GIVE SINGULARITY DIAGNOSTIC
      KORD1I = 5
      go to 2240
 1200 TSPECS(1) = TN + HH
      if (LEX .eq. 0) go to 1250
      if (HH * (TSPECS(1) - TMARKX) .lt. C0) go to 1250
      TSPECS(1) = TMARKX
      HH = TMARKX - TN
      LINC = 64
      if (LEX .gt. 0) go to 1240
      if ((LSC .lt. 4) .and. (HH / XI(1) .lt. CP3)) go to 1230
 1220 HH = CP875 * HH
      go to 1110
c     **** GIVE OUTPUT AT CURRENT TMARK (WITH EXTRAPOLATION)
 1230 KORD1I = -KMARK
      go to 1770
c     **** INTEGRATE TO TMARKX
 1240 KQICON = -1
c   TEST IF SUBROUTINE FOR COMPUTING INTEGRATION COEFF. SHOULD BE CALLED
 1250 continue
c++  Code for STIFF is inactive
c      IF ((KQMAXI .LT.KQICON) .OR. (KQMAXD.LT.KQDCON)) GO TO 1320
c++  Code for ~STIFF is active
      if (KQMAXI .lt. KQICON) go to 1320
c++  End
c   GO COMPUTE COEFFICIENTS REQUIRED FOR THE INTEGRATION
c     TEST IF STARTING
      if (LSC .lt. 7) go to 1310
 1260 KQMAXI = 2
c++  Code for STIFF is inactive
c      IF (METHOD) 1262,1270,1264
c 1262 KQMAXI=0
c 1264 KQMAXD=max(MAXDIF,2)
c      CALL SIVAHC
cc.  SET UP TO GO DO INITIALIZATION FOR CASE OF STIFF EQUATIONS
c      KORD1I=5
c      GO TO 1350
c++  End
c   INITIALIZE FOR EQUATIONS WHICH ARE NOT STIFF
 1270 KQMAXD = 0
      call SIVAHC
      J = NDTF
      do 1300 I = 1, NTE
c++  Code for STIFF is inactive
c         if (KORD(I + 3) .le. 0) go to 1290
c++  End
         KORD(I + 3) = 1
c     INITIALIZE THE DIFFERENCE TABLE
         if (LDT .eq. -4) F(J) = F(I)
         F(J + 1) = C0
         F(J + 2) = C0
 1290    continue
         J = J + NUMDT
 1300    continue
      if (LSC .eq. 5) go to 1340
      LSC = 7
      LDT = 1
      go to 1330
c   INTEGRATION IS NOT BEING STARTED
 1310 K = KORD(KEMAX + 3)
      SIGMAS = SIGMA(K)
      call SIVAHC
c     **** ADJUST EAVE
      TPS1 = BETA(K)
      if (TPS1 .gt. C1) TPS1 = CP5 * TPS1 + CP5
      EAVE = EAVE * TPS1 * (SIGMA(K) / SIGMAS)
c     TEST BELOW USED TO GET SAME RESULTS WITH/WITHOUT EXTRA EQUATIONS
      if (K .gt. KQICON) LSC = max(LSC, -3)
c END OF SPECIAL LOGIC FOR CASE WHEN INTEG. COEFF. ROUTINE IS CALLED
 1320 continue
c ********
c PREDICT Y
c ********
 1330 continue
c++  Code for ~ARGM is active
      call SIVAPR(Y, Y(NYNY), F, KORD)
c++  Code for ARGM is inactive
c      CALL SIVAPE
c++  End
c     GO GET PREDICTED DERIVATIVES
 1340 KORD1I = KPRED
c ********
c CALL SIVAF  (OR RETURN)
c ********
 1350 KORD2I = 0
 1360 KORD(1) = KORD1I
      KORD(2) = 0
      if (IOP13 .ne. 0) return
      call SIVAF(TSPECS(1), Y, F, KORD(1))
c     TEST FOR SPECIAL USER RETURN
 1380 if (KORD(1) .lt. 0) go to 2130
c     TEST FOR SPECIAL CASE
      if (KORD2I .ne. 0) go to 660
c ********
c TRANSFER CONTROL TO PROPER PLACE AFTER COMPUTING DERIVATIVES
c ********
      if (KORD1I - 2) 1400, 800, 1390
 1390 continue
c++  Code for VAREQ is active
      if (ICS - ICF) 1410, 1410, 760
c++  End
c ********
c PREPARE FOR CORRECTING, AND CORRECT Y
c ********
 1400 ITOLEP = 0
      ILGREP = 0
      IY = 1
      EIMAX = C1M5
      EMAX = C0
      KQMAXI = 2
      KQMAXS = 2
c++  Code for STIFF is inactive
c      IF (METHOD) 1404,1410,1406
c 1404 KQMAXI=0
c 1406 KQMAXD=2
c++  End
 1410 continue
c++  Code for ~ARGM is active
      call SIVACR(Y, F, KORD, F(NTOLF), KORD(IOP16))
c++  Code for ARGM is inactive
c      CALL SIVACE
c++  End
c     TEST IF ESTIMATED ERROR IS TOO BIG (OR IF DIAGNOSTIC CALLED FOR)
      if (EMAX .gt. EREP) if (EREP) 2210, 2210, 1670
 1420 continue
c++  Code for VAREQ is active
      if (IOP18 .ne. 0) go to 760
c++  End
 1430 KORD1I = 2
c     TEST IF NOISE APPEARS TO LIMIT PRECISION
      if (EMAX .lt. C0) go to 1470
c++  Code for ~STIFF is active
      if (LSC) 1450, 1360, 1610
c++  Code for STIFF is inactive
c      IF (LSC) 1450,1460,1610
c++  End
c     SET LSC=0 IF NOISE NO LONGER APPEARS TO LIMIT PRECISION
c     OR IF THE END OF THE STARTING PHASE HAS BEEN REACHED
 1450 LSC = 0
 1460 if (METHOD) 800, 1350, 1350
c ********
c NOISE APPEARS TO BE LIMITING PRECISION
c ********
 1470 continue
      if (LSC .le. 0) LSC = max(LSC - 1, -KQMAXS)
      if (LSC .eq. -1) go to 1460
      if (abs(EMAX) .lt. EXR) go to 1590
      LINC = -7
      TPS2 = (C1 + BETA(NOISEQ - 1)) ** NOISEQ
      if (SNOISE .lt. EEPS10 * TPS2) go to 1550
      TP = sign(EEPT75 * abs(TN) + OVTM75, HH)
      if (abs(TP) .gt. abs(HH)) go to 1550
      TSPECS(1) = TN + TP
      KORD1I = 0
c     **** GO TO BACK UP THE DIFFERENCES AND GET F(TSPECS(1))
      go to 1100
c     **** SOLUTION HAS BEEN INTERPOLATED AND F COMPUTED
 1490 continue
      KORD1I = 0
      LINC = LINC - 1
      if (LINC + 9) 1510, 1520, 1500
 1500 TSPECS(1) = TN + (TP + TP)
      TP1 = F(KEMAX)
      TP2 = F(NDTF + NUMDT * KEMAX - NUMDT)
      go to 1780
c     **** COMPUTE 2-ND DIFFERENCE AT CLOSE SPACED T VALUES
 1510 TP2 = TP3
 1520 TP3 = F(KEMAX)
      TPS1 = abs((TP3 - TP1) - (TP1 - TP2))
      if ((C16 * TPS1 * TPS2) .ge. DNOISE) if (LINC + 9) 1550, 870, 1550
 1530 continue
      TPS2 = CP25 * SNOISE / RBQ(NOISEQ)
      do 1540 K = 2, NUMDT
         TPS1 = TPS1 + TPS1
         RBQ(K) = max(TPS1, TPS2 * RBQ(K))
 1540    continue
      LINC = 0

cFTK Next two lines added 2009-10-15
      if (abs(EMAX) .lt. EREP) go to 1460
      LINC = -1

      HH = CP875 * HH
      go to 1040
c     **** SET UP TO GIVE NOISE DIAGNOSTIC
 1550 KORD1I = 6
      go to 2240
c     **** AFTER GIVING NOISE DIAGNOSTIC
 1560 KORD1I = 2
      if (KORD(2) .ge. 0) then
        TPS1 = EEPS10
cFTK Next line added 2009-10-15
        if (TPS1 .lt. .49E0 * RBQ(2)) go to 1530
      end if
c     **** SET NEW VALUE FOR OFFENDING TOL
      F(NTOLF + ITOLEP - 1) = FDAT(7)
      if (LINC + 7) 1180, 1570, 1180
 1570 LINC = 0
 1580 if (LSC) 1460, 1460, 1610
c     **** CHANGE HINCC AND ADJUST SIGMA( )
 1590 if (LSC .ne. -4) go to 1580
      if (HINCC .eq. C1P125) go to 1580
      TPS1 = C1P125 / HINCC
      TPS2 = 1.0E0
      do 1600 K = 2, IOP11
         TPS2 = TPS2 * TPS1
         SIGMA(K) = SIGMA(K) * TPS2
 1600    continue
      EAVE = EAVE * TPS1 ** (1-KORD(KEMAX+3))
      LINCD = 6
      LINCQ = 12
      HINCC = C1P125
      go to 1460
c   END OF CODE FOR CASE WHEN NOISE APPEARS TO LIMIT PRECISION
c ********
c SPECIAL LOGIC FOR STARTING THE INTEGRATION
c ********
 1610 if (LSC .eq. 1) go to 800
      LSC = LSC - 1
      if (LSC - 2) 1620, 1640, 1650
 1620 if (EIMAX .le. (CP0625*EAVE*(SIGMA(KQMAXS)/SIGMAS)*(BETA(KQMAXS+
     1   1))**2)) go to 800
 1630 KSSTRT = KSTEP + 2
c   TEST IF STEPSIZE IS TOO SMALL BEFORE ENDING STARTING PHASE
      if (abs(HH) .ge. HMIN) go to 1450
c     GIVE DIAGNOSTIC FOR STEPSIZE TOO SMALL AT END OF START
      KORD1I = 7
      go to 2240
c   SET LSC TO DO ONE DERIVATIVE EVAL. PER STEP
 1640 LSC = 1
      go to 800
c     TEST IF FIRST TIME THROUGH THE FIRST STEP
 1650 if (LSC .eq. 6) go to 1340
c     END STARTING PHASE IF CONVERGENCE OF CORRECTOR ITERATES TOO SLOW
      if (LDT .eq. -5) go to 1660
      LSC = min(LSC, 4)
      go to 800
 1660 LDT = 0
      if (LSC - 4) 1260, 1630, 1260
c END OF SPECIAL LOGIC FOR STARTING THE INTEGRATION
c ********
c ESTIMATED ERROR IS TOO BIG
c ********
 1670 if (BETA(2) - C1) 1690, 1730, 1680
 1680 HC = C1 / BETA(2)
      if (BETA(2) .ge. C1P125) go to 1740
 1690 if (BETA(2) .gt. CP3) go to 1730
c   REQUIRED STEPSIZE REDUCTION IS TOO RAPID -- GIVE A DIAGNOSTIC
      KORD1I = 4
      go to 2240
c
c     TEST KORD(2) AFTER ABOVE DIAGNOSTIC OR A DISCONTINUITY DIAGNOSTIC
 1700 continue
      if (KORD(2) .eq. 0) go to 1730
c  TEST IF SOLUTION MUST BE DUMPED BEFORE A RESTART
 1710 LINC = -1
c++  Code for DUMP is active
      if (IOP9 .eq. 0) go to 1750
      if (KIS .eq. 2) go to 1750
      LINC = -6
c.    GO DUMP SOLUTION BEFORE REPEATING THE STEP
 1720 KQMAXI = KQMXIS
c++  Code for DUMP & STIFF is inactive
c      KQMAXD=KQMXDS
c++  Code for DUMP is active
      call SIVABU(F, KORD)
      go to 900
c++  End
c   SET UP TO REPEAT THE STEP
 1730 HC = CP5
 1740 LINC = -1
      if (LSC .le. 3) go to 1030
c   RESTART THE INTEGRATION IF ERROR IS TOO BIG ON FIRST OR SECOND STEP
c LOOP TO SELECT A NEW INITIAL STEPSIZE
 1750 LSC = 7
 1755 HH = HH * CP5
      EMAX = EMAX * CP25
      if (EMAX .ge. CP3) go to 1755
      go to 1090
c   END OF SELECTING A NEW INITIAL STEPSIZE
c END OF LOGIC FOR CASE WHEN ESTIMATED ERROR IS TOO BIG
c ********
c INTEGRATION HAS REACHED AN OUTPUT POINT
c ********
 1760 if (KMARK .eq. 0) go to 1920
      KORD1I = min(KMARK, 5)
      KORD(3) = KMARK
      if (TSPECS(1) .eq. TMARK) go to 1790
 1770 TSPECS(1) = TMARK
 1780 call SIVAIN(TSPECS(1), Y, F, KORD)
 1790 continue
      if (KORD1I) 1800, 730, 1810
c   OUTPUT POINT IS OBTAINED BY EXTRAPOLATION
 1800 continue
c++  Code for EXTRAP is active
      KORD1I = -KORD1I
      KORD2I = -7
      KEXIT = 4
c.  TEST IF GSTOP-S ARE PRESENT
c++  Code for EXTRAP &  GSTOP is active
      if (NGTOT .eq. 0) go to 1820
      IGFLG = 4
      KEXIT = 2
      KORD1I = 7
      if (IOP7) 740, 1820, 740
c++  End
c ********
c CALL SIVAO  (OR RETURN)
c ********
 1810 KORD2I = 1
 1820 KORD(1) = KORD1I
      KORD(2) = 1
      if (IOP14 .ne. 0) return
      call SIVAO(TSPECS(1), Y, F, KORD(1))
c     TEST FOR SPECIAL USER RETURN OR FOR A RESTART
c++  Code for ~DUMP is inactive
c 1840 IF (KORD(1)) 2130,700,1880
c++  Code for DUMP is active
 1840 if (KORD(1) .gt. 0) go to 1880
 1850 if (IOP9 .eq. 0) go to 1870
c.    **** GO DUMP THE SOLUTION
      LINC = -7
      ITOLEP = KORD(1)
      IDAT(1) = KORD(2)
      NEPTOL = KORD1I
      if (LSC .ne. 8) go to 900
 1860 LINC = min(0, LINCD)
      KORD1I = NEPTOL
      KORD(1) = ITOLEP
      KORD(2) = IDAT(1)
 1870 if (KORD(1)) 2130, 700, 2100
c++  End
 1880 if (KORD2I .lt. 0) go to (2140, 1810, 1350, 2110, 720, 750, 680,
     1   710), -KORD2I
      if (KORD2I .eq. 0) go to 1380
c ********
c TRANSFER CONTROL TO PROPER PLACE AFTER OUTPUT
c ********
 1890 if (KORD1I - 5) 1910, 1930, 1900
 1900 if (KORD1I - 8) 840, 1130, 910
 1910 if (KORD1I - 3) 1920, 1930, 880
c   GET NEW TOUT
 1920 TOUT = TSPECS(1) + TSPECS(3)
c GET NEW TMARK (NEXT INDEP. VAR. OUTPUT POINT)
 1930 XP = TMARK
      K = KMARK
      TMARK = TOUT
      KMARK = 2
      LEX = 0
      if (IOP5) 1940, 1980, 1970
 1940 I = -IOP5
 1950 I = I + 3
      J1 = KORD(I - 3)
      if (J1) 1950, 1980, 1960
 1960 J2 = KORD(I - 2)
      L = KORD(I - 1)
      go to 1990
 1970 J1 = 5
      J2 = IOP5
      L = 0
      if (J2 .ge. J1) go to 1990
 1980 J1 = 4
      J2 = 4
      L = IOP3
c
c     **** LOOP TO SET NEW TMARK (AND TMARKX)
 1990 do 2060 J = J1, J2
c        **** TEST IF EXTRAPOLATION NOT POSSIBLE
         if (L .eq. 0) go to 2010
         LX = 2
         if (LEX) 2020, 2030, 2020
 2000    LEX = L
 2010    LX = 1
 2020    if (HH * (TSPECS(J) - TMARKA(LX))) 2030, 2060, 2060
 2030    if (J .eq. 4) go to 2050
         if (HH * (TSPECS(J) - XP)) 2060, 2040, 2050
 2040    if ((K .ge. J) .or. (K .eq. 3)) go to 2060
 2050    TMARKA(LX) = TSPECS(J)
         if (LX .eq. 2) go to 2000
         KMARK = J
 2060    continue
      if (IOP5 .lt. 0) go to 1950
      if (J1 .ne. 4) go to 1980
      if (KMARK .eq. 4) KMARK = 3
c     **** TEST IF NEW TMARK IS ACCEPTABLE
      if (HH * (XP - TMARK)) 2070, 2080, 2090
 2070 if (KORD2I - 1) 670, 840, 670
 2080 if (K .ne. KMARK) go to 2070
c++  Code for DUMP is active
      if (KORD1I .eq. 3) go to 1850
c++  Code for ~DUMP is inactive
c      IF (KORD1I .EQ. 3) GO TO 2100
c++  End
 2090 if (KORD1I .eq. 13) go to 2190
c SETUP TO INDICATE ERROR IN SPECIFICATION OF OUTPUT POINTS
      KORD1I = 2
      IDAT(2) = KMARK
      if (KMARK .le. 3) IDAT(2) = KMARK + 1
      FDAT(3) = TSPECS(IDAT(2))
      go to 2240
c     SET KORD1I=1 TO INDICATE THAT END OF INTEGRATION HAS BEEN REACHED
 2100 KORD1I = 1
c ********
c RETURN TO USER
c ********
 2110 KORD2I = -1
      KORD(1) = KORD1I
 2130 KORD(2) = -1
      return
c ********
c TRANSFER CONTROL TO PROPER PLACE AFTER RETURN TO USER
c ********
 2140 if (KORD1I - 2) 2150, 1560, 2160
 2150 KORD2I = 1
      go to 1930
 2160 if (KORD1I - 4) 1700, 2200, 2170
 2170 if (KORD1I - 13) 2180, 1930, 2190
 2180 if (abs(HH) .ge. HMIN) if (KORD1I - 11) 1030, 1450, 1030
      if (KORD(2) .eq. 0) if (KORD1I - 11) 1070, 800, 1070
c   ERROR MESSAGES HAVE BEEN IGNORED -- COMPUTATION CAN NOT CONTINUE
 2190 KORD1I = 1
      go to 2240
c
c        AFTER A DISCONTINUITY RETURN
 2200 LINC = -4
      if (KORD(2)) 1710, 1730, 1100
c ********
c PROBLEM ENCOUNTERED WHEN CORRECTING
c ********
 2210 if (LDIS .eq. 0) go to 2230
c           Extra checks when had a user specified discontinuity.
c++  Code for VAREQ is active
      if (IOP18 .ne. 0) go to 760
c++  End
 2220 KORD1I = 2
      LDIS = LDIS + 1
      TP = DISADJ / HH
      if (KIS .ge. 1000) then
         if (LDIS .eq. 2) then
            if (KQMAXS .le. 3) then
               LDIS = 0
               EREP = abs(EREP)
               TSPECS(2) = HH*min(min(TP, TP**2),
     1            (CP25 * EXR/EMAX)**.333333333E0)
               go to 720
            end if
            LINC = -5
            if (IOP9 .eq. 0) KIS = 1001
            go to 800
         end if
         if (IOP9 .eq. 0) KIS = KIS + 1
         if (KQMAXS .le. LDIS + 2) KIS = LDIS + 1
         LINC = min(LINC, LDIS-2)
      end if
      if (LDIS .gt. 2*KQMAXS) then
         EREP = abs(EREP)
         LDIS = 0
         if (EMAX .gt. EREP) go to 1670
         go to 1430
      end if
      if (TP .ge. HINCC**(LINC+2)) then
         if ((LDIS .ne. 3) .and. (TP .gt. real(KQMAXS))) LSC = 1
         EIMIN = CP5
         EAVE = EAVE * TP**8
      end if
      if (LSC .eq. 2) go to 1630
      if (EMAX .gt. EXR) go to 1730
      go to 1430
c
 2230 EREP = abs(EREP)
c++  Code for DUMP is active
      if (LINC .lt. -3) go to 1720
c++  End
c     BAD TOL
      KORD1I = 3
c ********
c ERROR PROCESSING
c ********
 2240 FDAT(1) = TN
      FDAT(2) = HH
      IDAT(1) = KSTEP
      ITOLEP = max(NEPTOL, -NEPTOL - 1)
      J = 3
      if (KORD1I .ge. 7) then
         J = 4
         FDAT(3) = HMIN
      end if
      if (KORD1I .le. 3) then
         if (KORD1I .lt. 3) then
            K = 8
         else
            MACT(9) = LTXTAL
            FDAT(3) = C0
            IDAT(2) = ITOLEP
            IDAT(3) = ITOLEP + NTOLF - 1
            K = 11
         end if
      else
         MACT(9) = LTXTAK
         FDAT(J) = EMAX
         IDAT(2) = KEMAX
         IDAT(3) = ITOLEP
         IDAT(4) = ITOLEP + NTOLF - 1
         FDAT(J+1) = F(IDAT(4))
         K = 14
         if (KORD1I .eq. 6) then
            K = 17
            IDAT(5) = IDAT(4)
            FDAT(7) = 32.E0 * abs(EMAX) * FDAT(J+1)
            FDAT(J+2) = FDAT(7)
         end if
         MACT(12) = LTXTAL
         if (NEPTOL .lt. 0) then
            MACT(12) = LTXTAM
            IDAT(6) = IDAT(4)
            IDAT(5) = IDAT(4) + 1
            FDAT(J+2) = F(IDAT(5))
            FDAT(J+3) = FDAT(J+1) * FDAT(J+2)
         end if
      end if
c Set the location for the first part of the message that varies, set
c the error severity, and the index number, print the error and
c return or stop.
      L = MLOC(KORD1I)
      MACT(6) = L / LOCM
      MACT(2) = (L - MACT(6) * LOCM) / 32
      KORD1I = mod(L, 32)
      MACT(3) = KORD1I
      MACT(K) = MERET
c--S Next line special: P=>S, X=>D
      call SMESS(MACT, MTXTAA, IDAT, FDAT)
      MACT(K) = MENTXT
      go to 2110
c
      end
c   End of SIVAA

      subroutine SIVABU(F, KORD)
c>> 1987-12-07 SIVABU Krogh   Initial code.
c
c THIS SUBROUTINE RESTORES THE DIFFERENCE TABLE TO ITS STATE
c AT THE BEGINNING OF THE CURRENT STEP.  IF THE INTEGRATION ORDER
c WAS INCREASED, IT IS REDUCED. THE COMMON ARRAY XI IS ALSO
c RESTORED TO ITS STATE AT THE BEGINNING OF THE STEP. IF THE
c STEPSIZE IS NOT BEING CHANGED, THE ARRAY V USED TO COMPUTE
c INTEGRATION COEFFICIENTS IS RESTORED.
c
      integer KORD(*)
      real             F(*)
c
      integer KDIM, MAXORD, MAXSTF
c++ Substitute for KDIM, MAXORD, MAXSTF below
      parameter (KDIM = 16, MAXORD = 2, MAXSTF = 1)
c--S Next line special: P=>D, X=>Q
      real             TN
      real             XI(KDIM)
c
c--S Next line special: P=>D, X=>Q
      real             TG(2), TGSTOP(2), TMARK, TMARKX, TOUT, TOLG
      real             ALPHA(KDIM), BETA(KDIM+1)
      real              D(MAXSTF+MAXORD,MAXORD), G(KDIM,MAXORD)
      real             V(KDIM+MAXORD)
      real             HC, HDEC, HINC, HINCC, HMAX, HMAXP9, HMIN
      real             FDAT(11)
c
      real             DS(MAXSTF+MAXORD, MAXORD), GS(KDIM)
      real             SIGMA(KDIM), RBQ(KDIM), DNOISE
      real             EAVE, EIMAX, EIMIN, EMAX, EREP, ROBND, SNOISE
c
      integer IOPST, KORDI, KQMAXD, KQMAXI, LDT, MAXDIF, MAXINT, NKDKO,
     1   NTE, NYNY, NDTF, NUMDT
      common / SIVASC / TN, XI, IOPST, KORDI, KQMAXD, KQMAXI, LDT,
     1   MAXDIF, MAXINT, NKDKO, NTE, NYNY, NDTF, NUMDT
c
      integer ICF,ICS,IGFLG,IGTYPE(2),IGSTOP(2),ILGREP,INGS,IOP3,IOP4,
     1   IOP5,IOP6,IOP7,IOP8,IOP9,IOP10,IOP11,IOP12,IOP13,IOP14,IOP15,
     2   IOP16,IOP17,IOP18,IOP19,IOP20,IOP21,IOP22,IOP21S,ITOLEP,IY,
     3   KEMAX,KIS,KMARK,KORD1I,KORD2I,KPRED,KQDCON,KQICON,KQMAXS,
     4   KQMXDS,KQMXIL,KQMXIP,KQMXIS,KSC,KSOUT,KSSTRT,KSTEP,LEX,LINC,
     5   LINCD,LINCQ,LSC,MAXKQD,MAXKQI,METHOD,NE,NEPTOL,NG,NGTOT,
     6   NOISEQ,NOUTKO,NTOLF,NY,IDAT(6)
      common /SIVAMC/ TG,TGSTOP,TMARK,TMARKX,TOUT,TOLG,HC,HDEC,HINC,
     1   HINCC,HMAX,HMAXP9,HMIN,ALPHA,BETA,D,G,V,DS,GS,SIGMA,RBQ,DNOISE,
     2   EAVE,EIMAX,EIMIN,EMAX,EREP,ROBND,SNOISE,FDAT,ICF,ICS,IGFLG,
     3   IGTYPE,IGSTOP,ILGREP,INGS,IOP3,IOP4,IOP5,IOP6,IOP7,IOP8,IOP9,
     4   IOP10,IOP11,IOP12,IOP13,IOP14,IOP15,IOP16,IOP17,IOP18,IOP19,
     5   IOP20,IOP21,IOP22,IOP21S,ITOLEP,IY,KEMAX,KIS,KMARK,KORD1I,
     6   KORD2I,KPRED,KQDCON,KQICON,KQMAXS,KQMXDS,KQMXIL,KQMXIP,KQMXIS,
     7   KSC,KSOUT,KSSTRT,KSTEP,LEX,LINC,LINCD,LINCQ,LSC,MAXKQD,MAXKQI,
     8   METHOD,NE,NEPTOL,NG,NGTOT,NOISEQ,NOUTKO,NTOLF,NY,IDAT
      save / SIVAMC / , / SIVASC /
c
      integer I, L, KQQ, J, K
      real             TPD, C0, C2
      parameter (C0 = 0.E0)
      parameter (C2 = 2.E0)
c ********* START OF EXECUTABLE CODE **********
c
c ********
c BEGIN LOOP TO BACK UP DIFFERENCE TABLES
c ********
      L = NDTF - 1
      do 2410 I = 1, NTE
         KQQ = KORD(I + 3)
c++  Code for STIFF is inactive
c         IF (KQQ) 2302,2400,2310
cc.           EQUATION IS STIFF
c 2302    IF (LINC.GE.0) GO TO 2310
c         IF (F(L+1+I)) 2306,2308,2304
cc.     ORDER WAS INCREASED, AND THUS MUST BE DECREASED (KQQ.LT.0)
c 2304    KQQ=KQQ+1
c         KORD(I+3) = KQQ
c         GO TO 2308
cc.     ORDER WAS DECREASED
c 2306    KQQ=KQQ-1
c 2308    KQQ=max(2,-KQQ)
c         GO TO 2350
c++  End
c     EQUATION IS NOT STIFF
 2310    if (KQQ .gt. 2) then
            if (F(L + KQQ) .eq. C0) then
c                 ORDER WAS INCREASED, AND THUS MUST BE DECREASED
               KQQ = KQQ - 1
               KORD(I + 3) = KQQ
            end if
         end if
         J = min(KQQ, KSC)
         KQMAXI = max(KQMAXI, KQQ)
         if (KQQ .ne. 1) F(L + KQQ + 1) = 0.E0
c           BACK UP FOR BACKWARD DIFFERENCES
         do 2360 K = 1, J
            F(L + K) = F(L + K) - F(L + K + 1)
 2360    continue
         if (KQQ .gt. KSC) then
c           BACK UP FOR MODIFIED DIVIDED DIFFERENCES
            do 2390 K = J+1, KQQ
               F(L + K) = (F(L+K) - F(L+K+1)) / BETA(K)
 2390       continue
         end if
 2400    F(L + KQQ + 1) = F(L + KQQ + 1) / BETA(KQQ + 1)
         L = L + NUMDT
 2410 continue
c END OF LOOP TO BACK UP DIFFERENCE TABLES
c ********
c BACK UP XI TO BEGINNING OF THE STEP
c ********
      I = KSC + 1
      if (I - IOP11 - 1) 2420, 2440, 2450
 2420 TPD = XI(1)
c                Check below needed when starting?
      if (TPD .eq. XI(2)) go to 2450
      do 2430 K = I, IOP11
 2430    XI(K - 1) = XI(K) - TPD
 2440 XI(IOP11) = C2 * XI(IOP11 - 1)
      if (IOP11 .ne. 2) XI(IOP11) = XI(IOP11) - XI(IOP11 - 2)
 2450 KQICON = -1
      ICF = NE
      ICS = 1
      LDT = 1
      return
      end
c   End of SIVABU

      subroutine SIVACO(ID, RD)
c>> 1987-12-07 SIVACO Krogh   Initial code.
c
c THIS SUBROUTINE RETURNS THE FOLLOWING DATA FROM COMMON
c ID(1) = KEMAX  =  INDEX OF EQUATION WITH LARGEST ERROR ESTIMATE
c ID(2) = KSTEP  =  CURRENT STEP NUMBER
c ID(3) = NUMDT  =  NUMBER OF DIFFERENCES USED FOR EACH EQUATION
c ID(4) =           RESERVED FOR FUTURE USE
c ID(5) =           RESERVED FOR FUTURE USE
c RD(1) = EMAX   =  MAX. RATIO OF ESTIMATED ERROR TO REQUESTED ERROR
c RD(2) =           RESERVED FOR FUTURE USE
c RD(3) =           RESERVED FOR FUTURE USE
c
      integer ID(5)
      real             RD(3)
c
      integer KDIM, MAXORD, MAXSTF
c++ Substitute for KDIM, MAXORD, MAXSTF below
      parameter (KDIM = 16, MAXORD = 2, MAXSTF = 1)
c--S Next line special: P=>D, X=>Q
      real             TN
      real             XI(KDIM)
c
c--S Next line special: P=>D, X=>Q
      real             TG(2), TGSTOP(2), TMARK, TMARKX, TOUT, TOLG
      real             ALPHA(KDIM), BETA(KDIM+1)
      real              D(MAXSTF+MAXORD,MAXORD), G(KDIM,MAXORD)
      real             V(KDIM+MAXORD)
      real             HC, HDEC, HINC, HINCC, HMAX, HMAXP9, HMIN
      real             FDAT(11)
c
      real             DS(MAXSTF+MAXORD, MAXORD), GS(KDIM)
      real             SIGMA(KDIM), RBQ(KDIM), DNOISE
      real             EAVE, EIMAX, EIMIN, EMAX, EREP, ROBND, SNOISE
c
      integer IOPST, KORDI, KQMAXD, KQMAXI, LDT, MAXDIF, MAXINT, NKDKO,
     1   NTE, NYNY, NDTF, NUMDT
      common / SIVASC / TN, XI, IOPST, KORDI, KQMAXD, KQMAXI, LDT,
     1   MAXDIF, MAXINT, NKDKO, NTE, NYNY, NDTF, NUMDT
c
      integer ICF,ICS,IGFLG,IGTYPE(2),IGSTOP(2),ILGREP,INGS,IOP3,IOP4,
     1   IOP5,IOP6,IOP7,IOP8,IOP9,IOP10,IOP11,IOP12,IOP13,IOP14,IOP15,
     2   IOP16,IOP17,IOP18,IOP19,IOP20,IOP21,IOP22,IOP21S,ITOLEP,IY,
     3   KEMAX,KIS,KMARK,KORD1I,KORD2I,KPRED,KQDCON,KQICON,KQMAXS,
     4   KQMXDS,KQMXIL,KQMXIP,KQMXIS,KSC,KSOUT,KSSTRT,KSTEP,LEX,LINC,
     5   LINCD,LINCQ,LSC,MAXKQD,MAXKQI,METHOD,NE,NEPTOL,NG,NGTOT,
     6   NOISEQ,NOUTKO,NTOLF,NY,IDAT(6)
      common /SIVAMC/ TG,TGSTOP,TMARK,TMARKX,TOUT,TOLG,HC,HDEC,HINC,
     1   HINCC,HMAX,HMAXP9,HMIN,ALPHA,BETA,D,G,V,DS,GS,SIGMA,RBQ,DNOISE,
     2   EAVE,EIMAX,EIMIN,EMAX,EREP,ROBND,SNOISE,FDAT,ICF,ICS,IGFLG,
     3   IGTYPE,IGSTOP,ILGREP,INGS,IOP3,IOP4,IOP5,IOP6,IOP7,IOP8,IOP9,
     4   IOP10,IOP11,IOP12,IOP13,IOP14,IOP15,IOP16,IOP17,IOP18,IOP19,
     5   IOP20,IOP21,IOP22,IOP21S,ITOLEP,IY,KEMAX,KIS,KMARK,KORD1I,
     6   KORD2I,KPRED,KQDCON,KQICON,KQMAXS,KQMXDS,KQMXIL,KQMXIP,KQMXIS,
     7   KSC,KSOUT,KSSTRT,KSTEP,LEX,LINC,LINCD,LINCQ,LSC,MAXKQD,MAXKQI,
     8   METHOD,NE,NEPTOL,NG,NGTOT,NOISEQ,NOUTKO,NTOLF,NY,IDAT
      save / SIVAMC / , / SIVASC /
c
      ID(1) = KEMAX
      ID(2) = KSTEP
      ID(3) = NUMDT
      RD(1) = EMAX
      return
      end
c   End of SIVACO

      subroutine SIVACR(Y, F, KORD, TOL, LGROUP)
c>> 1988-08-25 SIVACR Krogh   Fix bug in relative error test.
c>> 1988-01-15 SIVACR Krogh   Initial code.
c
c THIS SUBROUTINE
c   1. CORRECTS Y FOR EQUATIONS WHICH ARE NOT STIFF
c   2. ESTIMATES ERRORS
c   3. SELECTS INTEGRATION ORDERS
c   4. TESTS IF NOISE LIMITS THE PRECISION
c
c     Y = VECTOR OF PREDICTED VALUES ON ENTRY, AND OF CORRECTED
c         VALUES WHEN THE RETURN IS MADE.
c LGROUP= VECTOR INDICATING HOW ERROR TOLERANCES ARE TO BE GROUPED
c         (AND POSSIBLY HOW INTEGRATION ORDERS ARE TO BE GROUPED).
c   TOL = VECTOR CONTAINING ERROR TOLERANCES (AND POSSIBLY RELATIVE
c         ERROR FACTORS).
c     F = VECTOR GIVING PREDICTED DERIVATIVE VALUES AND DIFF. TABLES.
c    KD = VECTOR GIVING ORDERS OF THE DIFFERENTIAL EQUATIONS
c         (IF EQUATIONS HAVE DIFFERENT ORDERS).
c    KQ = VECTOR OF INTEGRATION ORDERS.
c
      integer LGROUP(*), KORD(*)
c--S Next line special: P=>D, X=>Q
      real             Y(*)
      real             TOL(*), F(*)
c
      integer KDIM, MAXORD, MAXSTF
c++ Substitute for KDIM, MAXORD, MAXSTF below
      parameter (KDIM = 16, MAXORD = 2, MAXSTF = 1)
c--S Next line special: P=>D, X=>Q
      real             TN
      real             XI(KDIM)
c
c--S Next line special: P=>D, X=>Q
      real             TG(2), TGSTOP(2), TMARK, TMARKX, TOUT, TOLG
      real             ALPHA(KDIM), BETA(KDIM+1)
      real              D(MAXSTF+MAXORD,MAXORD), G(KDIM,MAXORD)
      real             V(KDIM+MAXORD)
      real             HC, HDEC, HINC, HINCC, HMAX, HMAXP9, HMIN
      real             FDAT(11)
c
      real             DS(MAXSTF+MAXORD, MAXORD), GS(KDIM)
      real             SIGMA(KDIM), RBQ(KDIM), DNOISE
      real             EAVE, EIMAX, EIMIN, EMAX, EREP, ROBND, SNOISE
c
      integer IOPST, KORDI, KQMAXD, KQMAXI, LDT, MAXDIF, MAXINT, NKDKO,
     1   NTE, NYNY, NDTF, NUMDT
      common / SIVASC / TN, XI, IOPST, KORDI, KQMAXD, KQMAXI, LDT,
     1   MAXDIF, MAXINT, NKDKO, NTE, NYNY, NDTF, NUMDT
c
      integer ICF,ICS,IGFLG,IGTYPE(2),IGSTOP(2),ILGREP,INGS,IOP3,IOP4,
     1   IOP5,IOP6,IOP7,IOP8,IOP9,IOP10,IOP11,IOP12,IOP13,IOP14,IOP15,
     2   IOP16,IOP17,IOP18,IOP19,IOP20,IOP21,IOP22,IOP21S,ITOLEP,IY,
     3   KEMAX,KIS,KMARK,KORD1I,KORD2I,KPRED,KQDCON,KQICON,KQMAXS,
     4   KQMXDS,KQMXIL,KQMXIP,KQMXIS,KSC,KSOUT,KSSTRT,KSTEP,LEX,LINC,
     5   LINCD,LINCQ,LSC,MAXKQD,MAXKQI,METHOD,NE,NEPTOL,NG,NGTOT,
     6   NOISEQ,NOUTKO,NTOLF,NY,IDAT(6)
      common /SIVAMC/ TG,TGSTOP,TMARK,TMARKX,TOUT,TOLG,HC,HDEC,HINC,
     1   HINCC,HMAX,HMAXP9,HMIN,ALPHA,BETA,D,G,V,DS,GS,SIGMA,RBQ,DNOISE,
     2   EAVE,EIMAX,EIMIN,EMAX,EREP,ROBND,SNOISE,FDAT,ICF,ICS,IGFLG,
     3   IGTYPE,IGSTOP,ILGREP,INGS,IOP3,IOP4,IOP5,IOP6,IOP7,IOP8,IOP9,
     4   IOP10,IOP11,IOP12,IOP13,IOP14,IOP15,IOP16,IOP17,IOP18,IOP19,
     5   IOP20,IOP21,IOP22,IOP21S,ITOLEP,IY,KEMAX,KIS,KMARK,KORD1I,
     6   KORD2I,KPRED,KQDCON,KQICON,KQMAXS,KQMXDS,KQMXIL,KQMXIP,KQMXIS,
     7   KSC,KSOUT,KSSTRT,KSTEP,LEX,LINC,LINCD,LINCQ,LSC,MAXKQD,MAXKQI,
     8   METHOD,NE,NEPTOL,NG,NGTOT,NOISEQ,NOUTKO,NTOLF,NY,IDAT
      save / SIVAMC / , / SIVASC /
c.    SPECIFICATION OF ENVIRONMENTAL CONSTANTS.
      real             EEPS10, EEPS16, EROV10, EEPS2, EEPT75, EOVEP2
      real             OVTM75, OVD10
      common / SIVAEV / EEPS2, EEPT75, EOVEP2, OVTM75, OVD10, EEPS10,
     1   EEPS16, EROV10
      save / SIVAEV /
c
      integer L, I, KQL, KQN, KQD, JLGREP, J, K, ILGROR, ITOLOR, JLGROR,
     1   IORD, KOUTKO, KQLORD, LL, LKQMAX
      real             CM8, CM2, CMP5, C0, CQ3125, CP1, CP125, CP25, CP5
      real             CP75, CP8, CP9375, C1, C1P4, C2, C4, C10, C20
      real             C1000, C40
      parameter (CM8 = -8.E0)
      parameter (CM2 = -2.E0)
      parameter (CMP5 = -.5E0)
      parameter (C0 = 0.E0)
      parameter (CQ3125 = .03125E0)
      parameter (CP1 = .1E0)
      parameter (CP125 = .125E0)
      parameter (CP25 = .25E0)
      parameter (CP5 = .5E0)
      parameter (CP75 = .75E0)
      parameter (CP8 = .8E0)
      parameter (CP9375 = .9375E0)
      parameter (C1 = 1.E0)
      parameter (C1P4 = 1.4E0)
      parameter (C2 = 2.E0)
      parameter (C4 = 4.E0)
      parameter (C10 = 10.E0)
      parameter (C20 = 20.E0)
      parameter (C40 = 40.E0)
      parameter (C1000 = 1000.E0)
      real             TPP, HH, E, EI, EPS, ERCOEF, RND, RNOISE, S
      real             TP2, TPS1, TPS2, TPS3, TPS4, TPS5, TPS6, TPS7
      real             REF(4)
      real             EIBND(KDIM-1)
c++  Code for INTEGO is active
      real             TEMPA(4), TEMPAO(4)
c++  End
      save KOUTKO, LKQMAX
      equivalence (TPS1,TEMPA(1)), (TPS2,TEMPA(2)), (TPS3,TEMPA(3)),
     1   (TPS4, TEMPA(4))
      equivalence (G(1, 1), HH)
      integer MACT1(2), MACT2(12)
c             Parameters for Interface to MESS and SMESS
      integer MERET, METEXT, METABL
      parameter (MERET  =51)
      parameter (METEXT =53)
      parameter (METABL =55)
c ********* Error message text ***************
c[Last 2 letters of Param. name]  [Text generating message.]
cAA KSTEP=$(I6) T=$(E15.8) H=$(E12.5) LSC=$(I3) $C
c   EIMIN=$(E8.2) EAVE=$G KSC=$(I2) SIGMA($J)=$G $C
c   RQ=$(E11.5)$G$E
c   $
cAB I$HKQ$HLI$HE$HEI$HEPS$HF$H$H$H$H
c   HIGH ORDER PREDICTED DIFFERENCES$HRNOISE$HSTIFF$HBETA$E
      integer LTXTAA,LTXTAB
      parameter (LTXTAA=  1,LTXTAB=  1)
      character MTXTAA(1) * (104)
      character MTXTAB(1) * (88)
      data MTXTAA/'KSTEP=$(I6) T=$(E15.8) H=$(E12.5) LSC=$(I3) EIMIN=$(E
     *8.2) EAVE=$G KSC=$(I2) SIGMA($J)=$G RQ=$(E11.5)$G$E'/
      data MTXTAB/'I$HKQ$HLI$HE$HEI$HEPS$HF$H$H$H$HHIGH ORDER PREDICTED$
     * DIFFERENCES$HRNOISE$HSTIFF$HBETA$E'/
c
      data MACT1 / METEXT, MERET /
c (rr=repeat, t=3/5 for I/E format)  wwddtrr  wwddtrr  wwddtrr
      data MACT2 / METABL, 1, 0, 14, 0400201, 0300202, 0801503,
     1    1507501, 1103504, 0901501, 1002501, 1205501 /
c         wwddtrr  wwddtrr  wwddtrr  wwddtrr  wwddtrr
c          End of stuff for interface to message processor
c
      data REF(1), REF(2), REF(3), REF(4) / C1, CP9375, CP75, CP5 /
c++ Save data by elements if ~.C.
c++ Of next 20 lines, only the first KDIM-1 are active
      data EIBND(1) / .1E0 /
      data EIBND(2) / .1E0 /
      data EIBND(3) / .14E0 /
      data EIBND(4) / .19E0 /
      data EIBND(5) / .26E0 /
      data EIBND(6) / .36E0 /
      data EIBND(7) / .50E0 /
      data EIBND(8) / .69E0 /
      data EIBND(9) / .94E0 /
      data EIBND(10) / C1 /
      data EIBND(11) / C1 /
      data EIBND(12) / C1 /
      data EIBND(13) / C1 /
      data EIBND(14) / C1 /
      data EIBND(15) / C1 /
C     data EIBND(16) / C1 /
C     data EIBND(17) / C1 /
C     data EIBND(18) / C1 /
C     data EIBND(19) / C1 /
C     data EIBND(20) / C1 /
c
c++  Code for ARGM is inactive
c      RETURN
c      ENTRY SIVACE
c++  End
c ********
c START OF CODE
c ********
      L = NDTF - 1
      if (ICS .ne. 1) L = L + (ICS - 1) * NUMDT
      do 3340 I = ICS, ICF
         if (NKDKO .ne. 0) KORDI = KORD(NKDKO + I - 1)
         IY = IY + abs(KORDI)
         KQL = KORD(I + 3)
         KQN = abs(KQL)
         KQD = max(2, KQN)
c ********
c OBTAIN ERROR TOLERANCE SPECIFIED BY THE USER
c ********
         if (I .le. ILGREP) if (KQL) 2600, 3310, 2610
         ITOLEP = abs(ITOLEP) + 1
         EPS = TOL(ITOLEP)
         ILGREP = LGROUP(ITOLEP)
c   TEST IF SIMPLE ABSOLUTE ERROR TEST IS BEING USED
         if (ILGREP .gt. 0) go to 2580
         JLGREP = ILGREP
c     GET OLD RELATIVE ERROR FACTOR
         TPS6 = TOL(ITOLEP + 1)
         ILGREP = LGROUP(ITOLEP + 1)
         ITOLEP = -ITOLEP - 1
c
         if (JLGREP + 1) 2540, 2570, 2510
c   NO CHECK ON THE ERROR ESTIMATE IS TO BE MADE
 2510    if (EPS + C1) 2520, 2590, 2520
c   ERROR TOLERANCE IS SPECIFIED IMPROPERLY
 2520    KEMAX = I
         NEPTOL = ITOLEP
         LINC = -3
         EREP = -abs(EREP)
         return
c   COMPUTE NEW RELATIVE ERROR FACTOR
 2540    continue
         TPS1 = C0
         do 2550 J = I, ILGREP
            TPS1 = TPS1 + abs(F(J))
 2550       continue
         TPS1 = abs(HH) * TPS1 / real(ILGREP - I + 1)
         if (LSC .le. 2) go to 2560
c     ON FIRST 3 STEPS INCREASE TPS6 WHEN COMPUTING REL. ERROR FACTOR
         TPS6 = max(C4 * TPS1, TPS6)
c     ON 1-ST TIME THROUGH THE FIRST STEP, REL. ERR. FAC. IS NOT STORED
         if (LSC .eq. 7) go to 2570
 2560    continue
         TPS6 = max(TPS1, TPS6)
c   STORE NEW RELATIVE ERROR FACTOR
         TOL(-ITOLEP) = TPS6 * REF(-JLGREP - 1)
c   COMPUTE ABSOLUTE ERROR TOLERANCE
 2570    EPS = EPS * TPS6
 2580    if (EPS .le. C0) go to 2520
 2590    if (KQL) 2600, 3330, 2610
c END OF OBTAINING ERROR TOLERANCE
c ********
c OBTAIN INFORMATION USED FOR ERROR ESTIMATION, ORDER SELECTION, ETC.
c ********
c EQUATION IS STIFF
 2600    continue
c++  Code for STIFF is inactive
c      JS=abs(KORD(NJSKO+I-1))-1
c      JSI=JS
c      TPP=C0
c      TPS4=F(L+KQD+2)
c      TPS3=F(L+KQD+1)
c      TPS2=F(L+KQD)
c      TPS1=F(L+KQD-1)
c      IF (KQD.EQ.2) TPS1=Y(IY-1)
c      E=ABS(TPS3)+ABS(TPS4)
c      EI=E+ABS(TPS2)
c      RND=EI
c      IF (KORDI.GE.0) GO TO 2604
cc.    EQUATION IS IMPLICIT
c      JSI=JSI-1
c      IF (JSI.NE.0) GO TO 2604
c      IF (KORDI.EQ.-1) GO TO 2602
c      ERCOEF=GS(KQN+1)
c      GO TO 2606
c 2602 ERCOEF=.5E0*DS(KQD,1)
c      JSI=1
c      GO TO 2606
cc.    END OF SPECIAL CODE FOR IMPLICIT EQUATIONS
c 2604 ERCOEF = DS(KQD,JSI)
c 2606 ERCOEF = ABS(ERCOEF) / EPS
c      IF (LSC.LE.2)  GO TO 2710
c      IF (LSC-5) 2650,2710,2710
cc.  END OF CODE FOR STIFF EQUATIONS
c++  End
c
c EQUATION IS NOT STIFF
 2610    TPP = F(I) - F(L + 1)
         TPS3 = TPP
         TPS4 = TPP - F(L + KQD + 1)
         TPS2 = TPP + F(L + KQD)
         TPS1 = TPP + F(L + KQD - 1)
         E = abs(TPS3) + abs(TPS4)
         RND = E
         EI = E + abs(TPS2)
         ERCOEF = abs(GS(KQN + 1)) / EPS
         if (KQL .ge. 4) go to 2710
c   TEST IF STARTING OR IF INTEGRATION ORDER IS ONE
         if (LSC .le. 2) if (KQL - 2) 2660, 2710, 2710
c ********
c LOGIC ASSOCIATED WITH STARTING THE INTEGRATION
c ********
         TPS4 = C0
         if (LSC - 4) 2650, 2640, 2620
c FIRST STEP
 2620    E = E * CQ3125
         TPS3 = C0
         F(L + 4) = C0
         S = C0
c   TEST IF FIRST TIME THROUGH THE FIRST STEP
         if (LSC .eq. 7) go to 2690
c   COMPUTE S=ESTIMATE OF H * EIGENVALUE OF JACOBIAN = 2*(F(A)-F(B))/
c   (F(B)-F(C)) WHERE F(A)=CURRENT F(I), AND F(B) AND F(C) PRECEDING
c   VALUES OR ESTIMATES OF F(I)
         TPP = F(I) - F(L + 5)
         TPS4 = TPP
         E = C2 * abs(TPS4)
         if (S .ne. C0) S = (TPS4 + TPS4) / S
         if (S + CP125) 2630, 2700, 2700
c     SET LDT=-5  TO INDICATE POSSIBLE PROBLEMS DUE TO INSTABILITY
 2630    LDT = -5
         go to 2690
c   ADJUST CORRECTION MADE ON SECOND STEP
 2640    TPP = CP8 * TPP
c   ADJUST ESTIMATED ERRORS ON SECOND AND THIRD STEPS
 2650    E = abs(TPS3)
         RND = C4 * E
         go to 2710
c END OF SPECIAL LOGIC FOR STARTING
c ********
c INTEGRATION ORDER =1 IS TREATED AS A SPECIAL CASE
c ********
 2660    TPP = TPP + F(L + 2)
         if (BETA(2) .ge. C1P4) EI = EI * C1000
c   ESTIMATE NEW VALUE FOR S
         S = F(L + 4)
         if (S .eq. C0) go to 2680
         S = max(CM8, C2 * BETA(2) * (TPS1 - TPS2 - F(L + 5)) / S)
         if (S .ge. CMP5) go to 2670
c   MODIFY TPP (TO GET BETTER STABILITY CHARACTERISTICS)
         TPP = TPP * max(CP25, (CM2 - C2 * S) / (S * S))
 2670    TPS4 = TPS4 * abs(S)
 2680    E = CP25 * (E + abs(TPS4))
         EI = EI + abs(TPS4 * S)
c     STORE INFORMATION REQUIRED TO ESTIMATE S ON NEXT STEP
 2690    F(L + 4) = TPP
 2700    F(L + 5) = F(I)
c END OF SPECIAL CODE FOR INTEGRATION ORDER =1
c ********
c CODE FOR NOISE TEST AND GETTING ERROR ESTIMATE
c ********
 2710    E = E * ERCOEF
         RNOISE = C0
         if (EPS .lt. C0) go to 2810
         TPS5 = abs(F(L + 2)) + abs(F(I))
         if (TPS5 .eq. C0) go to 2760
 2720    RNOISE = RND / TPS5
         if (RNOISE .gt. RBQ(KQD)) if (RNOISE - C1) 2760, 2750, 2750
c   NOISE IS APPARENTLY SLOWING CONVERGENCE OF THE DIFFERENCES
c     REDUCE EI
         EI = RND
         TPS5 = abs(EEPS2 * Y(IY - 1)) / EPS
         if (TPS5 .lt. abs(E)) if (LSC) 2730, 2730, 2760
         E = TPS5
         RNOISE = C0
 2730    E = -abs(E)
         if (EIMIN .gt. CP1) EI = (C10 * EIMIN) * EI
c     COMPUTE REDUCTION TO BE MADE IN EI
         if (RNOISE .gt. (C20 * RBQ(KQD))) go to 2760
         K = -6 - LSC
 2740    if (K .le. 0) go to 2760
c     REDUCE EI WHEN NOISE APPARENTLY LIMITS PRECISION
         K = K - 1
         EI = CP5 * EI
         if (EI .gt. EIMIN) go to 2740
         go to 2760
 2750    TPS4 = 1.1E0 * RND
         TPS3 = RND
 2760    continue
c   TEST FOR STIFFNESS GOES HERE WHEN IMPLEMENTED
c *       INGREDIENTS OF TEST MAY INCLUDE --
c *       RNOISE, WHETHER (ABS(TPS4).GT.ABS(TPS3)),
c *       WHETHER EMAX IS INCREASING, RESULT OF TEST ON
c *       PREVIOUS STEPS, ETC.
c
c ********
c COMPUTE ERROR ESTIMATES AND INFORMATION FOR SELECTING THE STEPSIZE
c ********
         if (E .ge. abs(EMAX)) go to 2770
         if (-E .le. abs(EMAX)) go to 2780
         SNOISE = RNOISE
         DNOISE = RND
         NOISEQ = KQD
c   STORE PARAMETERS ASSOCIATED WITH LARGEST VALUE OF E
 2770    EMAX = E
         KEMAX = I
         NEPTOL = ITOLEP
c   DETERMINE HOW MUCH STEPSIZE CAN BE INCREASED
 2780    EI = EI * ERCOEF * SIGMA(KQD)
         EIMAX = max(EIMAX, EI)
         if (LINC .le. 0) go to 2810
         K = 0
 2790    if (EI .ge. min(EIMIN, EIBND(KQN))) go to 2800
         K = K + 1
         if (K .eq. LINC) go to 2810
         EI = EI * SIGMA(KQD)
         go to 2790
 2800    LINC = K
c END OF COMPUTING ERROR ESTIMATES
 2810    continue
c++  Code for ERRSTO is inactive
c      IF (IOP20 .EQ. 0) GO TO 780
cc.********
cc.STORE ERROR ESTIMATE (OPTIONAL)
cc.********
c      F(IOP20+I-1)=TPS3*GS(KQN+1)
cc.END OF STORING ERROR ESTIMATE
c++  Code for INTEGO | ERRSTO is active
         if (IOP19 .eq. 0) go to 3090
c.********
c.EQUATIONS ARE GROUPED TO USE SAME INTEGRATION METHOD (OPTIONAL)
c.********
c++  Code for INTEGO is active
         if (I .gt. 1) if (I - ILGROR) 2900, 2900, 2830
         ITOLOR = IOP19
 2830    JLGROR = KORD(ITOLOR)
         ITOLOR = ITOLOR + 1
         if (JLGROR .gt. 0) go to 2870
         ILGROR = KORD(ITOLOR)
         ITOLOR = ITOLOR + 1
         if (JLGROR + 1) 2840, 2850, 2890
 2840    if (JLGROR .lt. -2) if (KQD + JLGROR) 2850, 2880, 2880
c.INITIALIZE FOR ACCUMULATING VARIABLES USED IN ORDER SELECTION
 2850    IORD = I
         KQLORD = KQL
         do 2860 K = 1, 4
 2860       TEMPAO(K) = abs(TEMPA(K))
         go to 2930
c.ORDERS IN CURRENT GROUP CAN BE DIFFERENT
 2870    ILGROR = JLGROR
         go to 3090
c.ORDER IS NOT GOING TO BE CHANGED
 2880    JLGROR = 0
 2890    if (KQL) 3240, 3270, 3270
c.TAKE ACTION FOR EQUATION WHICH IS NOT THE FIRST IN THE GROUP
 2900    if (JLGROR) 2910, 2890, 3090
c.ACCUMULATE VARIABLES USED IN ORDER SELECTION
 2910    do 2920 K = 1, 4
 2920       TEMPAO(K) = TEMPAO(K) + abs(TEMPA(K))
c.    TEST IF THIS IS LAST EQUATION IN THE GROUP
 2930    if (I .ne. ILGROR) if (KQL) 3310, 3290, 3290
c.SET UP TO GO SELECT INTEGRATION ORDER
         KQL = 0
         do 2940 K = 1, 4
 2940       TEMPA(K) = TEMPAO(K)
         go to 3090
c.INTEGRATION ORDER HAS BEEN SELECTED
c++  Code for INTEGO | STIFF is active
 2950    continue
c++  Code for INTEGO is active
         KQL = KQLORD
         if (KQN - abs(KQL)) 2960, 2980, 3020
c.  TEST IF ORDER CAN BE DECREASED
 2960    if (JLGROR .ge. -2) if (KQL) 3010, 3040, 3040
c.    INTEGRATION ORDER WAS SELECTED OUTSIDE PERMITTED RANGE
 2970    KQN = abs(KQL)
c.    INTEGRATION ORDER IS NOT GOING TO BE CHANGED
 2980    if ((KQL .ne. 1) .or. (LSC .gt. 0)) if (KQL) 3030, 3040, 3040
c.    SET  4-TH ENTRY IN DIFFERENCE TABLES SO THAT STANDARD ADAMS
c.    METHOD IS USED WHEN KQL=1
 2990    do 3000 K = IORD, I
 3000       F(NDTF + K*NUMDT - NUMDT + 3) = C0
         go to 3270
c.  ORDER FOR STIFF EQUATION WAS REDUCED
 3010    continue
c++  Code for INTEGO & STIFF is inactive
c      IF (KQN.LT.JSI) GO TO 990
c      TPP=-C1
c      GO TO 1090
cc.  TEST IF ORDER CAN BE INCREASED
c++  Code for INTEGO is active
 3020    if (JLGROR .eq. -2) go to 2970
c++  Code for INTEGO & STIFF is inactive
c      IF (KQL.GE.0) GO TO 1140
c      IF ((JSI.NE.0).AND.(KQN.GT.(MAXKQD+JSI))) GO TO 990
c      TPP=C1
cc.  STORE RESULTS FOR STIFF EQUATIONS
c++  Code for INTEGO is active
 3030    continue
c++  Code for INTEGO & STIFF is inactive
c      DO 3035 K=IORD,I
c      KORD(K+3) = -KQN
c 3035 F(NDTF+K*NUMDT-NUMDT)=TPP
c      GO TO 3245
cc.  STORE RESULTS FOR EQUATIONS WHICH ARE NOT STIFF
c++  Code for INTEGO is active
 3040    LL = NDTF + NUMDT * IORD - NUMDT
         do 3080 J = IORD, I
            KORD(J + 3) = KQN
            if (KQN - KQL) 3050, 3070, 3060
 3050       F(LL + KQD - 1) = F(LL + KQD - 1) + (F(J) - F(LL))
            go to 3080
 3060       F(LL + KQN) = F(LL + KQD)
 3070       F(LL + KQD) = C0
 3080       LL = LL + NUMDT
         if (KQN - 1) 3270, 2990, 3270
c++  End
c.********
c.SELECT INTEGRATION ORDER
c.********
 3090    if (LSC .le. 0) go to 3120
c. SPECIAL ORDER SELECTION WHEN STARTING
         if (LSC - 3) 3110, 3210, 3100
 3100    if (LSC .eq. 5) if (S + .125E0) 3160, 3130, 3130
         if (LSC - 6) 3130, 3130, 3210
 3110    if (C40 * min(abs(TPS4), abs(TPS3)) .gt. abs(TPS2)) then
            if (EPS .ne. -C1) LSC = 2
         end if
         if (abs(TPS4) .lt. abs(TPS3)) if (C4 * abs(TPS4) - abs(TPS2))
     1      3130, 3130, 3210
c.  CHECK IF ORDER CAN BE INCREASED OR SHOULD BE DECREASED
 3120    TPS5 = ROBND * abs(TPS4)
         TPS6 = ROBND * (TPS5 + abs(TPS3))
         TPS7 = abs(TPS1) + abs(TPS2)
         if (TPS5 .ge. abs(TPS3)) go to 3140
         if (TPS6 .ge. TPS7) go to 3210
 3130    if (KQN .ge. MAXKQI) go to 3210
c.    INCREASE THE INTEGRATION ORDER
         KQN = KQN + 1
c++  Code for INTEGO | STIFF is active
         if (KQL) 3230, 2950, 3250
c++  Code for ~(INTEGO | STIFF) is inactive
c      GO TO 3250
c++  End
c.  CHECK IF ORDER SHOULD BE DECREASED
 3140    if (TPS6 .lt. TPS7) go to 3210
         if (TPS5 .lt. abs(TPS3 - TPS4)) go to 3210
         if ((TPS3.eq.TPS4) .and. (LSC.le.0)) go to 3210
         if (KQN - 2) 3210, 3160, 3180
 3160    KQN = 1
c++  Code for INTEGO | STIFF is active
         if (KQL) 3220, 2950, 3170
c++  End
c.    WHEN ORDER IS REDUCED TO 1 WITH ADAMS METHOD SET F(L+4)=0
 3170    F(L + 4) = C0
         go to 3260
c.    DECREASE THE INTEGRATION ORDER
 3180    KQN = KQN - 1
c++  Code for INTEGO | STIFF is active
         if (KQL) 3220, 2950, 3200
c++  End
 3200    F(L+KQD) = F(L+KQD) + TPP
         go to 3260
c   NO CHANGE IN INTEGRATION ORDER IS BEING MADE
 3210    continue
c++  Code for INTEGO is active
         if (KQL) 3240, 2950, 3270
c++  Code for ~INTEGO is inactive
c         TPS1 = EEPS10
c      GO TO 1530
c++  End
c END OF SELECTING INTEGRATION ORDER
c ********
c COMPUTE MAXIMUM INTEGRATION ORDERS AND SET NEW ONES (IF ANY)
c ********
c EQUATION IS STIFF
c     ORDER WAS DECREASED
c++  Code for INTEGO | STIFF is active
 3220    continue
c++  Code for STIFF is inactive
c      IF (KQN.LT.JSI) GO TO 3236
c      F(L+1)=-C1
c      GO TO 3233
cc.    ORDER WAS INCREASED
c++  Code for INTEGO |  STIFF  is active
 3230    continue
c++  Code for STIFF is inactive
c      IF ((JSI.NE.0).AND.(KQN.GT.(MAXKQD+JSI))) GO TO 3236
c      F(L+1)=C1
c 3233 KORD(I+3) = -KQN
c      GO TO 3245
c      ORDER WAS SET TO AN UNACCEPTABLE VALUE
c 3236 KQN=abs(KQL)
c      ORDER IS NOT BEING CHANGED
c++  Code for STIFF |  INTEGO is active
 3240    continue
c++  Code for STIFF is inactive
c      F(L+1)=C0
c 3245 IF (JSI.NE.0) KQMAXD=max(KQN,KQMAXD)
c      IF (JS.LT.abs(KORDI)) KQMAXI=max(KQN,KQMAXI)
c      GO TO 3290
c++  End
c EQUATION IS NOT STIFF
c     ORDER INCREASED
 3250    F(L + KQN + 1) = -F(L + KQD + 1)
         if (LSC .gt. 0) F(L + KQN + 1) = F(L + 1) - F(I)
c     ORDER CHANGED
 3260    KORD(I + 3) = KQN
 3270    KQMAXI = max(KQN, KQMAXI)
         if (EPS .gt. C0) KQMAXS = max(KQN, KQMAXS)
         F(L + KQD + 1) = C0
 3290    continue
         if (KQN .gt. KIS) go to 3310
c.********
c.DETERMINE IF TIME TO STORE SOLUTION (OPTIONAL)
c.********
         if (KIS .ge. 1000) then
            TP2 = max(1.5E0, real(KQN) * C2 ** (1001 - KIS)) * abs(TPS4)
 3295       if (TP2 .gt. abs(F(L+KQN))) then
               if (KQN .le. KQL) then
                  KQN = KQN - 1
                  if (KQN .gt. 1) go to 3295
                  KQN = 1
               end if
            end if
            KORD(I+3) = KQN
            if (I .eq. 1) LKQMAX = 0
            LKQMAX = max(KQN, LKQMAX)
            KQMAXI = LKQMAX
            if (KIS .eq. 1000) then
               if (I .eq. KEMAX) EMAX = real(8 + KQN**2) * abs(EMAX)
               go to 3325
            end if
c++  Code for DUMP is active
         else if ((E .ne. C0) .and. (EPS .gt. C0)) then
            if (IOP9 .gt. 0) if((abs(E)*real(KIS-KQN+2)**(KQN+1))-1.E-2)
     1         3310, 3310, 3300
 3300       KIS = -1
c++  End
         end if
 3310    continue
c ********
c CORRECT
c ********
         do 3320 K = 1, KORDI
c++  Code for ~{p,x} is active
            Y(IY - K) = Y(IY - K) + G(KQL + 1, K) * TPP
c++  Code for {p,x} is inactive
Cc--S Next line special: P=>D, X=>Q
C            Y(IY - K) = Y(IY - K) + real(G(KQL + 1, K)) * real(TPP)
c++  END
 3320    continue
c END OF CORRECTING
 3325 continue
c++  Code for OUTPUT is active
      if (IOP10 .gt. 0) then
         if (I .eq. 1) then
            IDAT(1) = KSTEP
            IDAT(2) = LSC
            IDAT(3) = KSC
            IDAT(4) = IOP11
            FDAT(1) = TN
            FDAT(2) = HH
            FDAT(3) = EIMIN
            FDAT(4) = EAVE
            FDAT(5) = SIGMA(IOP11)
            FDAT(6) = ROBND
            MACT2(3) = NTE
c--S Next line special: P=>S, X=>D
            call SMESS(MACT1, MTXTAA, IDAT, FDAT)
            KOUTKO = NOUTKO
         end if
         if (KOUTKO .ne. 0) then
            if (KORD(KOUTKO) .gt. 0) then
               if (I .lt. KORD(KOUTKO)) go to 3328
               KOUTKO = KOUTKO + 1
            else
               if (I .ge. abs(KORD(KOUTKO))) KOUTKO = KOUTKO + 1
            end if
         end if
         IDAT(1) = I
         IDAT(2) = KQL
         IDAT(3) = LINC
         FDAT(1) = E
         FDAT(2) = EI
         FDAT(3) = EPS
         FDAT(4) = F(I)
         FDAT(5) = TPS1
         FDAT(6) = TPS2
         FDAT(7) = TPS3
         FDAT(8) = TPS4
         FDAT(9) = RNOISE
         FDAT(10) = 0.E0
         if (KQL .eq. 1) FDAT(10) = S
         FDAT(11) = BETA(KQD)
c--S Next line special: P=>S, X=>D
         call SMESS(MACT2, MTXTAB, IDAT, FDAT)
 3328    if (I .eq. NTE) IOP10 = IOP10 - 1
      end if
c++  End
 3330    L = L + NUMDT
 3340    continue
      return
      end
c   End of SIVACR

      subroutine SIVAHC
c>> 1988-05-20 SIVAHC Krogh   Initial code.
c
c SUBROUTINE TO COMPUTE COEFFICIENTS REQUIRED FOR INTEGRATING
c ORDINARY DIFFERENTIAL EQUATIONS
c
      integer KDIM, MAXORD, MAXSTF
c++ Substitute for KDIM, MAXORD, MAXSTF below
      parameter (KDIM = 16, MAXORD = 2, MAXSTF = 1)
c--S Next line special: P=>D, X=>Q
      real             TN
      real             XI(KDIM)
c
c--S Next line special: P=>D, X=>Q
      real             TG(2), TGSTOP(2), TMARK, TMARKX, TOUT, TOLG
      real             ALPHA(KDIM), BETA(KDIM+1)
      real              D(MAXSTF+MAXORD,MAXORD), G(KDIM,MAXORD)
      real             V(KDIM+MAXORD)
      real             HC, HDEC, HINC, HINCC, HMAX, HMAXP9, HMIN
      real             FDAT(11)
c
      real             DS(MAXSTF+MAXORD, MAXORD), GS(KDIM)
      real             SIGMA(KDIM), RBQ(KDIM), DNOISE
      real             EAVE, EIMAX, EIMIN, EMAX, EREP, ROBND, SNOISE
c
      integer IOPST, KORDI, KQMAXD, KQMAXI, LDT, MAXDIF, MAXINT, NKDKO,
     1   NTE, NYNY, NDTF, NUMDT
      common / SIVASC / TN, XI, IOPST, KORDI, KQMAXD, KQMAXI, LDT,
     1   MAXDIF, MAXINT, NKDKO, NTE, NYNY, NDTF, NUMDT
c
      integer ICF,ICS,IGFLG,IGTYPE(2),IGSTOP(2),ILGREP,INGS,IOP3,IOP4,
     1   IOP5,IOP6,IOP7,IOP8,IOP9,IOP10,IOP11,IOP12,IOP13,IOP14,IOP15,
     2   IOP16,IOP17,IOP18,IOP19,IOP20,IOP21,IOP22,IOP21S,ITOLEP,IY,
     3   KEMAX,KIS,KMARK,KORD1I,KORD2I,KPRED,KQDCON,KQICON,KQMAXS,
     4   KQMXDS,KQMXIL,KQMXIP,KQMXIS,KSC,KSOUT,KSSTRT,KSTEP,LEX,LINC,
     5   LINCD,LINCQ,LSC,MAXKQD,MAXKQI,METHOD,NE,NEPTOL,NG,NGTOT,
     6   NOISEQ,NOUTKO,NTOLF,NY,IDAT(6)
      common /SIVAMC/ TG,TGSTOP,TMARK,TMARKX,TOUT,TOLG,HC,HDEC,HINC,
     1   HINCC,HMAX,HMAXP9,HMIN,ALPHA,BETA,D,G,V,DS,GS,SIGMA,RBQ,DNOISE,
     2   EAVE,EIMAX,EIMIN,EMAX,EREP,ROBND,SNOISE,FDAT,ICF,ICS,IGFLG,
     3   IGTYPE,IGSTOP,ILGREP,INGS,IOP3,IOP4,IOP5,IOP6,IOP7,IOP8,IOP9,
     4   IOP10,IOP11,IOP12,IOP13,IOP14,IOP15,IOP16,IOP17,IOP18,IOP19,
     5   IOP20,IOP21,IOP22,IOP21S,ITOLEP,IY,KEMAX,KIS,KMARK,KORD1I,
     6   KORD2I,KPRED,KQDCON,KQICON,KQMAXS,KQMXDS,KQMXIL,KQMXIP,KQMXIS,
     7   KSC,KSOUT,KSSTRT,KSTEP,LEX,LINC,LINCD,LINCQ,LSC,MAXKQD,MAXKQI,
     8   METHOD,NE,NEPTOL,NG,NGTOT,NOISEQ,NOUTKO,NTOLF,NY,IDAT
      save / SIVAMC / , / SIVASC /
c.    SPECIFICATION OF ENVIRONMENTAL CONSTANTS.
      real             EEPS10, EEPS16, EROV10, EEPS2, EEPT75, EOVEP2
      real             OVTM75, OVD10
      common / SIVAEV / EEPS2, EEPT75, EOVEP2, OVTM75, OVD10, EEPS10,
     1   EEPS16, EROV10
      save / SIVAEV /
c                 K - 1 + 1 / K  is equivalent to max(1, K-1)
      real             GG(MAXORD - 1 + 1/MAXORD), B(KDIM+MAXORD),
     1   W(KDIM+MAXORD)
      integer  K, N, J
      real             C0, CP1, CRBQI, CP5, CP5625, C1, C1P125
      parameter (C0 = 0.E0)
      parameter (CP1 = .1E0)
      parameter (CRBQI = .421875E0)
      parameter (CP5 = .5E0)
      parameter (CP5625 = .5625E0)
      parameter (C1 = 1.E0)
      parameter (C1P125 = 1.125E0)
c++  Code for STIFF is inactive
c      INTEGER          GODIF
c++  End
      real             TP1, TP2, HH, TEMP, TP
      equivalence (G(1, 1), HH)
c
      save GG, W
c
c  B(K)= 1/(K*(K+1))
c++ Save data by elements if ~.C.
c++ Of next 23 lines, only the first KDIM+MAXORD are active
      data B(1)  / 5.000000000000000000000000000000000000000E-1 /
      data B(2)  / 1.666666666666666666666666666666666666667E-1 /
      data B(3)  / 8.333333333333333333333333333333333333333E-2 /
      data B(4)  / 5.000000000000000000000000000000000000000E-2 /
      data B(5)  / 3.333333333333333333333333333333333333333E-2 /
      data B(6)  / 2.380952380952380952380952380952380952381E-2 /
      data B(7)  / 1.785714285714285714285714285714285714286E-2 /
      data B(8)  / 1.388888888888888888888888888888888888889E-2 /
      data B(9)  / 1.111111111111111111111111111111111111111E-2 /
      data B(10) / 9.090909090909090909090909090909090909091E-3 /
      data B(11) / 7.575757575757575757575757575757575757576E-3 /
      data B(12) / 6.410256410256410256410256410256410256410E-3 /
      data B(13) / 5.494505494505494505494505494505494505495E-3 /
      data B(14) / 4.761904761904761904761904761904761904762E-3 /
      data B(15) / 4.166666666666666666666666666666666666667E-3 /
      data B(16) / 3.676470588235294117647058823529411764706E-3 /
      data B(17) / 3.267973856209150326797385620915032679739E-3 /
      data B(18) / 2.923976608187134502923976608187134502924E-3 /
C     data B(19) / 2.631578947368421052631578947368421052632E-3 /
C     data B(20) / 2.380952380952380952380952380952380952381E-3 /
C     data B(21) / 2.164502164502164502164502164502164502165E-3 /
C     data B(22) / 1.976284584980237154150197628458498023715E-3 /
C     data B(23) / 1.811594202898550724637681159420289855072E-3 /
c
c ********
c START OF CODE
c ********
c     SET STEP NUMBER OF METHOD
c++  Code for STIFF is inactive
c      IOP11 = MIN(max(KQMAXI,KQMAXD) + 1), KDIM)
c++  Code for ~STIFF is active
      IOP11 = MIN(KQMAXI + 1, KDIM)
c++  End
c     TEST IF STEPSIZE WAS CHANGED
      if (KQICON .ge. 0) go to 3510
c ********
c STEPSIZE JUST CHANGED
c ********
c     SET CONSTANTS DEPENDING ON NEW STEPSIZE
      KQMXIL = KQMAXI
      TP1 = HH
      GG(1) = TP1 * TP1
      G(1, 2) = GG(1) * CP5
      if (MAXINT .le. 2) go to 3450
      do 3440 K = 3, MAXINT
         GG(K - 1) = G(1, K - 1) * TP1
         G(1, K) = GG(K - 1) / real(K)
 3440    continue
c     SET CONSTANTS INDICATING STEP CHANGE
 3450 KQICON = 0
c++  Code for STIFF is inactive
c      KQDCON=0
c++  End
      KQMXIP = 1
      KSC = 1
      if (LSC .lt. 7) go to 3490
c     SPECIAL SET-UP OF CONSTANTS ON THE VERY FIRST STEP
      HINCC = C1P125
      LINCD = 6
      LINCQ = 12
      if (HINC .gt. C0) go to 3460
      LINCD = -2
      LINC = -2
      ROBND = C1
 3460 SIGMA(1) = 1.0E0
      BETA(1) = C1
      do 3470 N = 1, IOP11
c++  Code for STIFF is inactive
c      D(1,N)=C0
c++  End
         XI(N) = TP1
         ALPHA(N) = C1
         BETA(N + 1) = C1
         SIGMA(N + 1) = real(N + 1) * SIGMA(N) * HINCC
 3470    continue
      TEMP = EEPS16
      RBQ(1) = C1
      RBQ(2) = CP1
      TP = CRBQI
c     **** IN THE LOOP BELOW RBQ(K) IS COMPUTED TO BE
c          APPROXIMATELY (3/4 ** ((K-1) ** 2 - 1) / 10
c          .5625 = (3/4) ** 2    TP = (3/4) ** (2*K -3)
      do 3480 K = 3, KDIM
         TEMP = TEMP + TEMP
         RBQ(K) = max(TEMP, RBQ(K - 1) * TP)
 3480    TP = TP * CP5625
      go to 3560
c     SET-UP AFTER THE FIRST STEP
 3490 TP2 = XI(1)
      XI(1) = TP1
      BETA(2) = TP1 / TP2
      K = 2
      if (HINCC .eq. HINC) go to 3540
      if ((LSC .ne. 0) .or. ((KSTEP-KSSTRT-KQMAXS) .lt. 10)) go to 3540
      HINCC = C1
      LINCD = 0
 3500 LINCD = LINCD + 1
      HINCC = HINCC * HINC
      if (HINCC .lt. 2.E0) go to 3500
      LINC = (LINC * (LINCD + LINCD)) / LINCQ
      LINCQ = LINCD + LINCD
      HINCC = HINC
      go to 3540
c END OF LOGIC FOR CASE WHEN STEPSIZE JUST CHANGED
c     TEST IF MAXIMUM INTEGRATION ORDER DID NOT INCREASE
 3510 if (KQMAXI .gt. KQMXIL) then
c ********
c INTEGRATION ORDER WAS INCREASED -- GET NEW V'S
c ********
         KQMXIL = KQMAXI
         KQMXIP = KQMXIL + MAXINT
         K = KQMXIP
         V(K) = B(K)
         if (KQICON .eq. 1) go to 3530
c     if (KQICON .eq. K) KQICON = KQICON - 1 --- Removed 1999-08-19
         do 3520 N = 2, KQICON
            K = K - 1
 3520       V(K) = V(K) - ALPHA(N) * V(K + 1)
c END OF GETTING NEW V'S
      else
         IOP11 = max(IOP11, KQMXIL+1)
      end if
 3530 if (IOP11 .le. KSC) go to 3560
c ********
c COMPUTE PARAMETERS WHICH ARE STILL CHANGING AS A RESULT OF
c A CHANGE IN THE STEPSIZE
c ********
      TP2 = XI(KSC)
c     UPDATE CONSTANT STEP COUNTER
      KSC = KSC + 1
      K = KSC
      BETA(K) = C1
 3540 continue
      TEMP = HINCC
c
c   LOOP TO COMPUTE NEW VALUES OF PARAMETERS
      do 3550 N = K, IOP11
         TP1 = TP2 + HH
         TP2 = XI(N)
         XI(N) = TP1
         ALPHA(N) = HH / TP1
         BETA(N + 1) = BETA(N) * (TP1 / TP2)
         TEMP = max(TEMP, real(N) * (ALPHA(N) * HINCC))
         SIGMA(N) = SIGMA(N - 1) * TEMP
 3550    continue
      if (IOP11 .ne. KDIM) XI(IOP11 + 1) = TP2 + HH
c END OF CODE FOR COMPUTING PARAMETERS WHICH ARE STILL CHANGING
c
 3560 if (KQICON .ge. KQMXIP) go to 3690
c ********
c COMPUTE INTEGRATION COEFFICIENTS WHICH ARE STILL CHANGING
c ********
      KQMXIL = max(KQMAXI, KQMXIL)
      KQMXIP = KQMXIL + MAXINT
      J = KQMXIP - KQICON
      N = KQICON + 1
      KQICON = N
      if (N .ne. 1) go to 3580
c INITIALIZE V AND W
      do 3570 K = 1, J
         V(K) = B(K)
 3570    W(K) = V(K)
      go to 3600
c UPDATE V AND INITIALIZE W
 3580 if (N .eq. KDIM) go to 3690
      do 3590 K = 1, J
         V(K) = V(K) - ALPHA(N) * V(K + 1)
 3590    W(K) = V(K)
c SET TRANSFER FOR LOOP BELOW DEPENDING ON VALUE OF MAXINT
 3600 continue
      go to 3660
c
 3640 J = J - 1
c INNER LOOP FOR COMPUTING INTEGRATION COEFFICIENTS
      do 3650 K = 1, J
 3650    W(K) = W(K) - ALPHA(N) * W(K + 1)
c     STORE INTEGRATION COEFFICIENTS
 3660 G(N + 1, 1) = HH * W(1)
      GS(N + 1) = G(N + 1, 1) - G(N, 1)
c++  Code for MAXORD >= 2 is active
      if (MAXINT .ge. 2) then
         G(N + 1, 2) = GG(1) * W(2)
c++  Code for MAXORD >= 3 is inactive
c        if (MAXINT .gt. 2) then
c           DO 3665 K=3,MAXINT
c3665          G(N+1,K)=GG(K-1)*W(K)
c        end if
c++  Code for MAXORD >= 2 is active
      end if
c++  End
      N = N + 1
      if (N .le. KQMXIL) go to 3640
c END OF COMPUTING INTEGRATION COEFFICIENTS
c
 3690 continue
c++  Code for STIFF is inactive
c      IF (KQDCON.GT.KQMAXD) GO TO 4662
cc.********
cc.COMPUTE DIFFERENTIATION COEFFICIENTS WHICH ARE STILL CHANGING
cc.********
cc.SET TRANSFER FOR LOOP BELOW, DEPENDING ON VALUE OF MAXDIF
c++  Code for STIFF & MAXORD >= 2 is inactive
c      IF (MAXDIF-2) 3693,3692,3691
c 3691 ASSIGN 3696 TO GODIF
c      GO TO 3694
c 3692 ASSIGN 3698 TO GODIF
c      GO TO 3694
c 3693 ASSIGN 3699 TO GODIF
c++  Code for STIFF is inactive
c 3694 KQDCON=KQDCON+1
cc.LOOP FOR COMPUTING DIFFERENTIATION COEFFICIENTS
c      DO 3699 N=KQDCON,KQMAXD
c      DS(N+1,2)=C1/XI(N)
c      D(N+1,1)=DS(N+1,2)+D(N,1)
c      DS(N+1,1)=DS(N+1,2)/D(N+1,1)
c++  Code for STIFF & MAXORD >= 2 is inactive
c      GO TO GODIF, (3696,3698,3699)
c 3696 CONTINUE
c++  Code for STIFF & MAXORD >= 3 is inactive
c      DO 3697 K=3,MAXDIF
c      DS(N+1,K)=D(N,K-2) * (K-1)/XI(N)
c 3697 D(N+1,K-1)=DS(N+1,K) + D(N,K-1)
c++  Code for STIFF is inactive
c 3698 CONTINUE
c++  Code for STIFF & MAXORD >= 2 is inactive
c      D(N+1,MAXDIF)=D(N,MAXDIF) + D(N,MAXDIF-1) * (MAXDIF)/XI(N)
c++  Code for STIFF is inactive
c 3699 CONTINUE
c++  End
c
c END OF COMPUTING DIFFERENTIATION COEFFICIENTS
      return
      end
c   End of SIVAHC

      subroutine SIVAIN(T, Y, F, KORD)
c>> 1988-01-14 SIVAIN Krogh   Initial code.
c
c  SUBROUTINE TO DO INTERPOLATION FOR VARIABLE ORDER INTEG. ROUTINE
c
      integer KORD(*)
c--S Next line special: P=>D, X=>Q
      real             T(*), Y(*)
      real             F(*)
      integer KDIM, MAXORD
c++ Substitute for KDIM, MAXORD below
      parameter (KDIM = 16, MAXORD = 2)
c--S Next line special: P=>D, X=>Q
      real             TN
      real             XI(KDIM)
c
      integer IOPST, KORDI, KQMAXD, KQMAXI, LDT, MAXDIF, MAXINT, NKDKO,
     1   NTE, NYNY, NDTF, NUMDT
      common / SIVASC / TN, XI, IOPST, KORDI, KQMAXD, KQMAXI, LDT,
     1   MAXDIF, MAXINT, NKDKO, NTE, NYNY, NDTF, NUMDT
      save / SIVASC /
      integer I, ICI, IDT, INTERP, INTEG, INTEGZ, IY, IYI, IYN, IYNI, J,
     1    K, KQMXI, KQMXS, KQQ, L, N
      real             C0, C1, C2
      parameter (C0 = 0.E0)
      parameter (C1 = 1.E0)
      parameter (C2 = 2.E0)
      real             C(KDIM+MAXORD-1), ETA(KDIM)
      real             GAMMA(KDIM)
      real             TP1, HI
      real             CSUM(KDIM+MAXORD-1)
c--S Next line special: P=>D, X=>Q
      real             XP1
      logical LNOTM1
c
c              Stuff for processing error messages
      integer IDAT(1)
      real             FDAT(6)
      integer MENTXT, MERET, MEEMES, METEXT
      parameter (MENTXT =23)
      parameter (MERET  =51)
      parameter (MEEMES =52)
      parameter (METEXT =53)
      integer MACT(8)
c ********* Error message text ***************
c[Last 2 letters of Param. name]  [Text generating message.]
cAA SIVAIN$B
cAB Interpolating at T(1)=$F with $B
cAC TN=$F, T(2)=$F and H=$F.  T(1) must be in [$F, $F].$E
cAD internal variable LDT = $I.  Interpolation not allowed now.$E
      integer LTXTAA,LTXTAB,LTXTAC,LTXTAD
      parameter (LTXTAA=  1,LTXTAB=  9,LTXTAC= 41,LTXTAD= 94)
      character MTXTAA(1) * (154)
      data MTXTAA/'SIVAIN$BInterpolating at T(1)=$F with $BTN=$F, T(2)=$
     *F and H=$F.  T(1) must be in [$F, $F].$Einternal variable LDT = $I
     *.  Interpolation not allowed now.$E'/
c
c                      1 2 3 4       5 6       7      8
      data MACT / MEEMES,0,0,0, MENTXT,0, METEXT, MERET /
c
c++  Code for ARGM is inactive
c      ENTRY SIVAIE
c++  End
c ********
c START OF CODE -- CHECK ON STATE OF DIFFERENCE TABLE
c ********
      L = LDT
      if (L) 3710, 3730, 3780
 3710 if (L + 2) 4170, 3730, 3720
 3720 if (MAXINT .ge. 0) L = 1
      go to 3840
c ********
c UPDATE DIFFERENCE TABLE TO START OF NEXT STEP
c ********
 3730 K = NDTF
      do 3770 I = 1, NTE
         KQQ = KORD(I + 3)
         if (KQQ .le. 0) go to 3760
c EQUATION IS NOT STIFF
         TP1 = F(I) - F(K)
c LOOP TO DO UPDATING
         N = K + max(abs(KQQ), 2)
         do 3750 J = K, N
 3750       F(J) = F(J) + TP1
 3760    continue
         K = K + NUMDT
 3770    continue
      LDT = 1
      if (L .ne. 0) return
c END OF UPDATING DIFFERENCE TABLE
c ********
c INITIALIZE FOR COMPUTATION OF COEFFICIENTS
c ********
 3780 INTERP = 0
      HI = T(1) - TN
      GAMMA(1) = HI / XI(1)
      if (GAMMA(1)) 3790, 3800, 3810
 3790 if (GAMMA(1) .ge. -C1) go to 3820
      INTERP = 1
      if (abs(HI) - abs(T(2))) 3820, 3820, 4180
 3800 INTERP = 2 - KQMAXI
      go to 3820
 3810 if (GAMMA(1) .gt. C2) if (LDT - 2) 4180, 3820, 3820
 3820 KQMXI = KQMAXI + INTERP - 1
c++  Code for STIFF is inactive
c      KQMXS=max(KQMXI,KQMAXD)
c++  Code for ~STIFF is active
      KQMXS = KQMXI
c++  End
      do 3830 N = 2, KQMXS
 3830    GAMMA(N) = (HI + XI(N-1)) / XI(N)
 3840 LNOTM1 = L .ne. -1
      INTEG = MAXINT
      if (INTEG .le. 0) if (INTEG + MAXDIF) 4160, 3950, 3950
c ********
c COMPUTE INTEGRATION COEFFICIENTS
c ********
c     INITIAL SET-UP
c         COMPUTE INITIAL C VALUES
      do 3850 N = 1, INTEG
         C(N) = HI / real(N)
 3850 continue
      I = INTEG + 1
      INTEG = INTEG + KQMXI
      do 3860 N = I, INTEG
         C(N) = C(N - 1) * (real(N - MAXINT) / real(N))
 3860    continue
c         COMPUTE ETA'S
      do 3870 N = 1, KQMXI
 3870    ETA(N) = HI / XI(N)
c         COMPUTE C(K)'S TO CORRESPOND TO G(K-MAXINT+1,MAXINT),
c         K=MAXINT, MAXINT+1,..., MAXINT+KQMXI-1
      I = INTEG
 3880 J = INTEG
      INTEG = J - 1
      if (INTEG .le. MAXINT) go to 3900
      do 3890 N = J, I
 3890    C(N) = ETA(N - INTEG) * C(N) + C(N - 1)
      go to 3880
 3900 do 3910 N = J, I
 3910    C(N) = ETA(N - INTEG) * C(N)
c         END OF COMPUTING  G(---,MAXINT)
      INTEGZ = 0
      go to 3940
c         COMPUTE C(K)-S TO CORRESPOND TO G(K-INTEG+1,INTEG),
c         K=INTEG+1,INTEG+2,..., INTEG+KQMXI
 3920 do 3930 N = 1, KQMXI
 3930    C(INTEG+N) = GAMMA(N)*C(INTEG+N-1) - ETA(N)*C(INTEG+N)
 3940 ICI = INTEG - 1
      go to 4020
c END OF COMPUTING INTEGRATION COEFFICIENTS
c ********
c COMPUTE COEFFICIENTS FOR INTERPOLATION
c ********
 3950 C(1) = C1
      ICI = 0
      do 3960 N = 1, KQMXS
 3960    C(N + 1) = GAMMA(N) * C(N)
      if (INTEG + 1) 3970, 3990, 4010
c END OF COMPUTING INTERPOLATION COEFFICIENTS
c
c     SET-UP TO COMPUTE DIFFERENTIATION COEFFICIENTS REQUIRED
c     IN ORDER TO GET COEFFICIENTS ACTUALLY USED
 3970 INTEG = 0
      ICI = 1
 3980 INTEG = INTEG - 1
      if (INTEG .eq. MAXINT) ICI = 0
c ********
c COMPUTE DIFFERENTIATION COEFFICIENTS
c ********
 3990 INTERP = max(INTERP, 0)
      TP1 = real(-INTEG)
      C(1) = TP1 * C(1) / XI(-INTEG)
      J = KQMAXD + INTEG
      do 4000 N = 1, J
 4000    C(N + 1) = (TP1*C(N)) / XI(N - INTEG) + GAMMA(N - INTEG) * C(N)
c     C(N) NOW CORRESPONDS TO THE DIFFERENTIAL COEFFICIENT
c          D(N-INTEG,-INTEG)
 4010 INTEGZ = INTEG
      if (ICI .ne. 0) go to 3980
c END OF COMPUTING DIFFERENTIATION COEFFICIENTS
c ********
c BEGINNING OF LOOP TO DO
c         INTEGRATION       (INTEG.GT.0)
c         INTERPOLATION     (INTEG.EQ.0)
c         DIFFERENTIATION   (INTEG.LT.0)
c TO THE POINT INDICATED BY T.
c ********
c     SET UP INITIAL INDICES
 4020 if (NYNY .lt. 0) then
         IY = -NYNY
         IYNI = NYNY + ICI + 1
         if (LDT .eq. 2) then
            CSUM(ICI+1) = C(ICI+1)
            do 4025 J = ICI+2, INTEG+KQMXI
            CSUM(J) = CSUM(J-1) + C(J)
 4025       continue
         end if
      else
         IY = 1
         IYNI = NYNY + ICI - 1
      end if
      IDT = NDTF - INTEGZ
      do 4140 I = 1, NTE
         if (NKDKO .ne. 0) KORDI = KORD(NKDKO + I - 1)
         IY = IY + abs(KORDI)
         KQQ = KORD(I + 3)
c         GET INDEX OF HIGHEST ORDER DIFFERENCE TO BE USED
         K = max(abs(KQQ) + INTERP, 2)
         IYI = -INTEG
         if (KQQ) 4030, 4130, 4040
c EQUATION IS STIFF
 4030    continue
c++  Code for STIFF is inactive
c      JS=abs(KORD(NJSKO+I-1))-1
c      IYI=IYI-JS
c      IF(LNOTM1) IF (IYI) 4034,4032,4130
c      IF (KORDI.LT.0) IYI=IYI+1
c      IYI=IYI+MAXINT-abs(KORDI)
c      IF (IYI) 4034,4130,4130
cc.      IF EQUATION IS IMPLICIT DO NOT COMPUTE AN F
c 4032 IF (KORDI.LT.0) GO TO 4130
cc.      TEST IF INTEG TOO BIG FOR THIS EQUATION
c 4034 IF (abs(KORDI).LT.-IYI) GO TO 4130
c      IYI=IYI+IY
c      IYN=IYI+IYNI
cc. COMPUTE INNER PRODUCT FOR STIFF EQUATIONS
c      IF (INTEGZ.EQ.0) GO TO ???
cc.    DIFFERENTIATING
c      TP1 = C0
c      DO 4036 J = K+INTEGZ, 1, -1
c         TP1 = TP1 + C(J) * F(IDT+J-1)
c 4036 CONTINUE
cc.    TEST WHETHER TO STORE RESULT IN Y OR F
c      IF (IYI-IY) 4080, 4090, 4080
cc.    INTEGRATING OR INTERPOLATING
c      TP1 = C0
c      DO 4037 J = ICI + K, ICI + 2, -1
c         TP1 = TP1 + C(J) * F(IDT+J-ICI-1)
c 4037 CONTINUE
c      IF (INTEG.EQ.0) GO TO 4120
c      TP1=TP1 + C(ICI+1)*Y(IYN+1)
c++  End
         go to 4100
c END OF SPECIAL CODE FOR STIFF EQUATIONS
c
c EQUATION IS NOT STIFF
 4040    if (LNOTM1) if (IYI) 4050, 4060, 4130
         IYI = IYI + MAXINT - KORDI
         if (IYI .ge. 0) go to 4130
c       TEST IF INTEG TOO BIG FOR THIS EQUATION
 4050    if (KORDI .lt. -IYI) go to 4130
 4060    IYI = IYI + IY
         IYN = IYI + IYNI
c  COMPUTE INNER PRODUCT FOR EQUATION WHICH IS NOT STIFF
         XP1 = C0
         if (LDT .eq. 2) then
            if (KQQ .ne. KQMAXI) XP1 = CSUM(K+INTEGZ+ICI) *
     1         F(IDT+INTEGZ+NUMDT-1)
         end if
         do 4070 J = K + INTEGZ + ICI, ICI + 1, -1
            XP1 = XP1 + C(J) * F(IDT - ICI - 1 + J)
 4070       continue
         if (INTEG) 4080, 4090, 4100
c STORE FINAL RESULT IN Y WHEN DIFFERENTIATING
 4080    continue
         Y(IYI) = XP1
         go to 4130
c STORE INTERPOLATED VALUE IN F (OR STIFF DIFFERENTIATION)
 4090    F(I) = XP1
         go to 4130
c PICK UP EXTRA STUFF TO ADD TO INNER PRODUCT WHEN INTEGRATING
 4100    K = ICI
         if (K .eq. 0) go to 4120
 4110    continue
         XP1 = C(K) * (XP1 + Y(IYN))
         IYN = IYN - 1
         K = K - 1
         if (K .ne. 0) go to 4110
c STORE FINAL RESULT IN Y WHEN INTEGRATING (OR STIFF INTERPOLATION)
 4120    Y(IYI) = XP1 + Y(IYN)
 4130    continue
         IDT = IDT + NUMDT
 4140    continue
c
      INTEG = INTEG - 1
      if (INTEG .ge. -MAXDIF) if (INTEG) 3990, 3950, 3920
 4160 return
c ********
c ERROR PROCESSING
c ********
 4170 MACT(2) = 68
      MACT(3) = 11
      MACT(6) = LTXTAD
      IDAT(1) = LDT
      go to 4190
 4180 MACT(2) = 28
      MACT(3) = 1
      MACT(6) = LTXTAC
      FDAT(2) = TN
      FDAT(3) = T(2)
      FDAT(4) = XI(1)
      FDAT(5) = TN - T(2)
      FDAT(6) = TN + C2 * XI(1)
      if (XI(1) .lt. 0) then
         FDAT(5) = FDAT(6)
         FDAT(6) = TN - T(2)
      end if
 4190 FDAT(1) = T(1)
c--S Next line special: P=>S, X=>D
      call SMESS(MACT, MTXTAA, IDAT, FDAT)
      if (MACT(2) .lt. 50) go to 3820
      return
      end
c   End of SIVAIN

      subroutine SIVAOP(IOPT, FOPT)
c>> 1987-12-07 SIVAOP Krogh   Initial code.
c
c  SUBROUTINE TO SET UP OPTIONS FOR DIFFERENTIAL EQUATION  PACKAGE -IVA
      real             FOPT(*)
      integer IOPT(*)
c
      integer KDIM, MAXORD, MAXSTF
c++ Substitute for KDIM, MAXORD, MAXSTF below
      parameter (KDIM = 16, MAXORD = 2, MAXSTF = 1)
c--S Next line special: P=>D, X=>Q
      real             TN
      real             XI(KDIM)
c
c--S Next line special: P=>D, X=>Q
      real             TG(2), TGSTOP(2), TMARK, TMARKX, TOUT, TOLG
      real             ALPHA(KDIM), BETA(KDIM+1)
      real              D(MAXSTF+MAXORD,MAXORD), G(KDIM,MAXORD)
      real             V(KDIM+MAXORD)
      real             HC, HDEC, HINC, HINCC, HMAX, HMAXP9, HMIN
      real             FDAT(11)
c
      real             DS(MAXSTF+MAXORD, MAXORD), GS(KDIM)
      real             SIGMA(KDIM), RBQ(KDIM), DNOISE
      real             EAVE, EIMAX, EIMIN, EMAX, EREP, ROBND, SNOISE
c
      integer IOPST, KORDI, KQMAXD, KQMAXI, LDT, MAXDIF, MAXINT, NKDKO,
     1   NTE, NYNY, NDTF, NUMDT
      common / SIVASC / TN, XI, IOPST, KORDI, KQMAXD, KQMAXI, LDT,
     1   MAXDIF, MAXINT, NKDKO, NTE, NYNY, NDTF, NUMDT
c
      integer ICF,ICS,IGFLG,IGTYPE(2),IGSTOP(2),ILGREP,INGS,IOP3,IOP4,
     1   IOP5,IOP6,IOP7,IOP8,IOP9,IOP10,IOP11,IOP12,IOP13,IOP14,IOP15,
     2   IOP16,IOP17,IOP18,IOP19,IOP20,IOP21,IOP22,IOP21S,ITOLEP,IY,
     3   KEMAX,KIS,KMARK,KORD1I,KORD2I,KPRED,KQDCON,KQICON,KQMAXS,
     4   KQMXDS,KQMXIL,KQMXIP,KQMXIS,KSC,KSOUT,KSSTRT,KSTEP,LEX,LINC,
     5   LINCD,LINCQ,LSC,MAXKQD,MAXKQI,METHOD,NE,NEPTOL,NG,NGTOT,
     6   NOISEQ,NOUTKO,NTOLF,NY,IDAT(6)
      common /SIVAMC/ TG,TGSTOP,TMARK,TMARKX,TOUT,TOLG,HC,HDEC,HINC,
     1   HINCC,HMAX,HMAXP9,HMIN,ALPHA,BETA,D,G,V,DS,GS,SIGMA,RBQ,DNOISE,
     2   EAVE,EIMAX,EIMIN,EMAX,EREP,ROBND,SNOISE,FDAT,ICF,ICS,IGFLG,
     3   IGTYPE,IGSTOP,ILGREP,INGS,IOP3,IOP4,IOP5,IOP6,IOP7,IOP8,IOP9,
     4   IOP10,IOP11,IOP12,IOP13,IOP14,IOP15,IOP16,IOP17,IOP18,IOP19,
     5   IOP20,IOP21,IOP22,IOP21S,ITOLEP,IY,KEMAX,KIS,KMARK,KORD1I,
     6   KORD2I,KPRED,KQDCON,KQICON,KQMAXS,KQMXDS,KQMXIL,KQMXIP,KQMXIS,
     7   KSC,KSOUT,KSSTRT,KSTEP,LEX,LINC,LINCD,LINCQ,LSC,MAXKQD,MAXKQI,
     8   METHOD,NE,NEPTOL,NG,NGTOT,NOISEQ,NOUTKO,NTOLF,NY,IDAT
      save / SIVAMC / , / SIVASC /
c.    SPECIFICATION OF ENVIRONMENTAL CONSTANTS.
      real             EEPS10, EEPS16, EROV10, EEPS2, EEPT75, EOVEP2
      real             OVTM75, OVD10
      common / SIVAEV / EEPS2, EEPT75, EOVEP2, OVTM75, OVD10, EEPS10,
     1   EEPS16, EROV10
      save / SIVAEV /
c
      integer IOPTS(23), INCOP(22), IOPTC(23), I, IA, J, K, LIOPT, MULTJ
      real             CMP75, C0, CP25, CP3, CP5, CP625, CP75, CP875,
     1   CP9, C1, C1P125, C2, C4, C10, C16
      parameter (CMP75 = (-.75E0))
      parameter (C0 = 0.E0)
      parameter (CP25 = .25E0)
      parameter (CP3 = .3E0)
      parameter (CP5 = .5E0)
      parameter (CP625 = .625E0)
      parameter (CP75 = .75E0)
      parameter (CP875 = .875E0)
      parameter (CP9 = .9E0)
      parameter (C1 = 1.E0)
      parameter (C1P125 = 1.125E0)
      parameter (C2 = 2.E0)
      parameter (C4 = 4.E0)
      parameter (C10 = 10.E0)
      parameter (C16 = 16.E0)
      external R1MACH
      real             R1MACH
      equivalence (IOPTC(3), IOP3)
      save IOPTS, LIOPT
c
c                      Declarations for error message processing.
c
      integer MECONT, MERET, MEEMES, MEIVEC
      parameter (MECONT =50)
      parameter (MERET  =51)
      parameter (MEEMES =52)
      parameter (MEIVEC =57)
      integer MACT(7), MACT1(5)
c
c ********* Error message text ***************
c[Last 2 letters of Param. name]  [Text generating message.]
cAA SIVAOP$B
cAB Error in IOPT() specifications: IOPT =$E
cAC HMIN = $F is > HMAX = $F.$E
      integer LTXTAA,LTXTAB,LTXTAC
      parameter (LTXTAA= 1,LTXTAB= 9,LTXTAC=49)
      character MTXTAA(1) * (76)
      data MTXTAA/'SIVAOP$BError in IOPT() specifications: IOPT =$EHMIN$
     * = $F is > HMAX = $F.$E'/
c **** End of text generated by pmess
c                      1   2  3   4       5  6      7
      data MACT / MEEMES, 88, 24, 0, MEIVEC, 0, MERET /
      data MACT1 / MEEMES, 28, 24, LTXTAC, MERET /
c
c                      IOP4       IOP17
      data IOPTS / 3*0, 500000, 12*0, 1, 6*0 /
c
c                  1  2    3  8  9 10 11 12  13  16   17 21 22
      data INCOP / 1, 3, 5*2, 1, 2, 3, 1, 2, 3*1, 3, 4*2, 2, 2 /
c
c ********* START OF EXECUTABLE CODE ***********************
c
      MULTJ = 1
      K = 1
 4200 I = IOPT(K)
      IA = abs(I)
c 1 and 6 lines below need 21 changed if more options are added.
      if (IA .le. 21) if (I) 4220, 4520, 4280
      if (IA .ne. 1111) go to 4490
      if (I .lt. 0) then
        MULTJ = -1
        K = K + 1
        go to 4200
      end if
      IOPT(2) = LIOPT
c
c     ****  INITIALIZE FOR STARTING A NEW INTEGRATION
      do 4210 J = 3, 23
 4210    IOPTC(J) = IOPTS(J)
      KSOUT = IOPTS(4)
      KMARK = 1 - IOPTS(1)
      KORDI = IOPTS(17)
      NKDKO = max(-KORDI, 0)
      IOPST = IOPTS(22)
      go to 4260
c
c     **** SET A NOMINAL VALUE
 4220 IOPTS(IA) = 0
      if (IA .eq. 12) go to 4420
      if (IA - 2) 4400, 4240, 4230
 4230 if (IA .eq. 4) IOPTS(4) = 500000
      if (IA .eq. 21) TOLG = 0.E0
      go to 4390
c
c     **** SET ALL OPTIONS TO THEIR NOMINAL VALUES
 4240 IA = 1
      IOPTS(1) = 0
      do 4250 J = 3, 22
         IOPTS(J) = 0
         IOPTC(J) = 0
 4250 continue
      IOPTS(4) = 500000
      IOPTC(4) = IOPTS(4)
      IOPTS(17) = 1
      TOLG = 0.E0
 4260 NGTOT = IOPTS(7) + max(IOPTS(6), 0)
      if (IOPTS(12) .eq. 0) go to 4420
 4270 return
c
c     **** SET SPECIFIED OPTION
 4280 J = IOPT(K + 1)
      if (INCOP(IA) - 2) 4290, 4330, 4300
c     **** OPTION INVOLVES NO EXTRA PARAMETERS
 4290 IOPTS(IA) = 1
      if (IA - 2) 4400, 4400, 4390
c     **** TAKE CARE OF SECOND EXTRA PARAMETER
 4300 if (IA .ne. 10) go to 4310
      NOUTKO = IOPT(K + 2)
      if (NOUTKO) 4500, 4350, 4350
 4310 if (IA .ne. 16) go to 4320
      NTOLF = IOPT(K + 2)
      if (NTOLF) 4500, 4500, 4350
 4320 if (J .eq. 3) then
        if (KMARK .ne. 3) then
           if (XI(1)*(FOPT(IOPT(K+2)) - TMARK) .ge. C0) go to 4400
        end if
      end if
      TMARK = FOPT(IOPT(K+2))
      KMARK = J
      go to 4400
c     **** TAKE CARE OF FIRST EXTRA PARAMETER
 4330 continue
      if (IA .eq. 12) go to 4410
      if (IA .eq. 4) KSOUT = J
      if (IA .eq. 21) TOLG = FOPT(J)
 4350 IOPTS(IA) = J * MULTJ
      if (abs(IA - 7) .gt. 1) go to 4360
c     **** SET SPECIAL PARAMETERS FOR GSTOP-S
      IGFLG = 0
      NGTOT = IOPTS(7) + max(IOPTS(6), 0)
c     **** TEST FOR ERROR
      if (J .gt. 500) go to 4500
 4360 if (J .gt. 0) go to 4390
      if ((IA .eq. 5) .or. (IA .eq. 17)) go to 4390
      if (J + 1) 4500, 4370, 4380
 4370 if (IA .eq. 7) go to 4500
 4380 if ((IA .eq. 4) .or. (IA .eq. 11) .or. (IA .ge. 16)) go to 4500
c     **** STORE SAVED VALUE IN COMMON
 4390 IOPTC(IA) = IOPTS(IA)
c
c     **** INCREMENT K TO GET NEXT OPTION
 4400 K = K + INCOP(IA)
      go to 4200
c
c ******* SET UP INFORMATION FOR CHANGING STEPSIZE *********
c
c     **** TEST IF VALUES ARE ALREADY SET
 4410 if (IOPTS(12) .ne. 0) go to 4430
c     **** SET NOMINAL VALUES FOR VARIABLES ONLY SET ONCE
 4420 EREP = CP3
c     **** SET NOMINAL VALUES FOR STEPSIZE CONTROL AND ENV. CONSTANTS
      EEPS2 = R1MACH(4)
      EEPS16 = C16 * EEPS2
      EEPS10 = CP625 * EEPS16
      EEPT75 = EEPS2 ** CP75
      EEPS2 = EEPS2 + EEPS2
      OVD10 = R1MACH(2)
      EROV10 = C10 / OVD10
      EOVEP2 = OVD10 * EEPS2
      OVTM75 = OVD10 ** CMP75
      OVD10 = OVD10 / C10
      HINC = C2
      HDEC = CP5
      HMIN = EROV10
      HMAX = OVD10
      if (I .ne. 12) go to 4470
 4430 IOPTS(12) = J * MULTJ
      if (J) 4450, 4470, 4460
c     **** SET UP TO GIVE USER COMPLETE STEPSIZE CONTROL
 4450 EREP = C1 / EROV10
      HINC = -C2
      IOP8 = 1
c## Recent code 12/16/94
      LINCD = -2
      LINC = -2
      ROBND = C1
c## Recent code 9/6/2001
      TOLG = 0.E0
c## End of recent code
      go to 4480
c     **** SET USER VALUES FOR STEPSIZE CONTROL
 4460 if (FOPT(J) .ne. C0) HINC = max(C1P125, min(FOPT(J),C4))
      if (FOPT(J + 1) .ne. C0) HDEC = min(CP875, max(FOPT(J + 1), CP25))
      if (FOPT(J + 2) .ne. C0) HMIN = FOPT(J + 2)
      if (FOPT(J + 3) .ne. C0) HMAX = FOPT(J + 3)
      if ((HMIN .gt. HMAX) .or. (HMAX .le. 0.E0)) then
         call SMESS(MACT1, MTXTAA, IOPT, FOPT(J+2))
         KORD2I = -4
      end if
 4470 KQICON = -1
 4480 HMAXP9 = HMAX * CP9
      if (I - 1111) 4400, 4270, 4400
c
c ***************** ERROR  IN  IOPT ************************
c
 4490 IA = 1
 4500 MACT(6) = K + INCOP(IA) - 1
      call MESS(MACT, MTXTAA, IOPT)
      KORD1I = 24
      KORD2I = -4
c Usual return with no error is here.
 4520 LIOPT = K
      return
      end
c   End of SIVAOP

      subroutine SIVAPR(Y, YN, F, KORD)
c>> 1988-01-13 SIVAPR Krogh   Initial code.
c
c THIS SUBROUTINE
c   1. UPDATES THE DIFFERENCE TABLE FROM THE PREVIOUS STEP (IF NOT
c      DONE ALREADY).
c   2. PREDICTS WHAT THE VALUES OF THE DEPENDENT VARIABLES, Y, AND
c      THE DIFFERENCE TABLE, DT, WILL BE AT THE END OF THE CURRENT STEP.
c
c   Y = VECTOR OF PREDICTED VALUES COMPUTED BY THIS SUBROUTINE.
c   YN= VECTOR OF VALUES OF Y COMPUTED ON THE LAST STEP.
c   F = VECTOR OF DERIVATIVE VALUES.
c   DT= ARRAY CONTAINING DIFFERENCE TABLES.
c   KD= VECTOR GIVING ORDERS OF THE DIFFERENTIAL EQUATIONS (IF
c       EQUATIONS HAVE DIFFERENT ORDERS).
c   KQ= VECTOR OF INTEGRATION ORDERS.
c
      integer KORD(*)
c--S Next line special: P=>D, X=>Q
      real             Y(*), YN(*)
      real             F(*)
c
      integer KDIM, MAXORD, MAXSTF
c++ Substitute for KDIM, MAXORD, MAXSTF below
      parameter (KDIM = 16, MAXORD = 2, MAXSTF = 1)
c--S Next line special: P=>D, X=>Q
      real             TN
      real             XI(KDIM)
c
c--S Next line special: P=>D, X=>Q
      real             TG(2), TGSTOP(2), TMARK, TMARKX, TOUT, TOLG
      real             ALPHA(KDIM), BETA(KDIM+1)
      real              D(MAXSTF+MAXORD,MAXORD), G(KDIM,MAXORD)
      real             V(KDIM+MAXORD)
      real             HC, HDEC, HINC, HINCC, HMAX, HMAXP9, HMIN
      real             FDAT(11)
c
      real             DS(MAXSTF+MAXORD, MAXORD), GS(KDIM)
      real             SIGMA(KDIM), RBQ(KDIM), DNOISE
      real             EAVE, EIMAX, EIMIN, EMAX, EREP, ROBND, SNOISE
c
c.    SPECIFICATION OF ENVIRONMENTAL CONSTANTS.
      real             EEPS10, EEPS16, EROV10, EEPS2
      real             EEPT75, EOVEP2, OVTM75, OVD10
      common / SIVAEV / EEPS2, EEPT75, EOVEP2, OVTM75, OVD10, EEPS10,
     1   EEPS16, EROV10
      save / SIVAEV /
      integer IOPST, KORDI, KQMAXD, KQMAXI, LDT, MAXDIF, MAXINT, NKDKO,
     1   NTE, NYNY, NDTF, NUMDT
      common / SIVASC / TN, XI, IOPST, KORDI, KQMAXD, KQMAXI, LDT,
     1   MAXDIF, MAXINT, NKDKO, NTE, NYNY, NDTF, NUMDT
c
      integer ICF,ICS,IGFLG,IGTYPE(2),IGSTOP(2),ILGREP,INGS,IOP3,IOP4,
     1   IOP5,IOP6,IOP7,IOP8,IOP9,IOP10,IOP11,IOP12,IOP13,IOP14,IOP15,
     2   IOP16,IOP17,IOP18,IOP19,IOP20,IOP21,IOP22,IOP21S,ITOLEP,IY,
     3   KEMAX,KIS,KMARK,KORD1I,KORD2I,KPRED,KQDCON,KQICON,KQMAXS,
     4   KQMXDS,KQMXIL,KQMXIP,KQMXIS,KSC,KSOUT,KSSTRT,KSTEP,LEX,LINC,
     5   LINCD,LINCQ,LSC,MAXKQD,MAXKQI,METHOD,NE,NEPTOL,NG,NGTOT,
     6   NOISEQ,NOUTKO,NTOLF,NY,IDAT(6)
      common /SIVAMC/ TG,TGSTOP,TMARK,TMARKX,TOUT,TOLG,HC,HDEC,HINC,
     1   HINCC,HMAX,HMAXP9,HMIN,ALPHA,BETA,D,G,V,DS,GS,SIGMA,RBQ,DNOISE,
     2   EAVE,EIMAX,EIMIN,EMAX,EREP,ROBND,SNOISE,FDAT,ICF,ICS,IGFLG,
     3   IGTYPE,IGSTOP,ILGREP,INGS,IOP3,IOP4,IOP5,IOP6,IOP7,IOP8,IOP9,
     4   IOP10,IOP11,IOP12,IOP13,IOP14,IOP15,IOP16,IOP17,IOP18,IOP19,
     5   IOP20,IOP21,IOP22,IOP21S,ITOLEP,IY,KEMAX,KIS,KMARK,KORD1I,
     6   KORD2I,KPRED,KQDCON,KQICON,KQMAXS,KQMXDS,KQMXIL,KQMXIP,KQMXIS,
     7   KSC,KSOUT,KSSTRT,KSTEP,LEX,LINC,LINCD,LINCQ,LSC,MAXKQD,MAXKQI,
     8   METHOD,NE,NEPTOL,NG,NGTOT,NOISEQ,NOUTKO,NTOLF,NY,IDAT
      save / SIVAMC / , / SIVASC /
      real             C0
      parameter (C0 = 0.E0)
c
      integer I, INTEG, INTEGS, J, K, KQQ, L, N
      real             TEMP(KDIM)
      real             TP1
c--S Next line special: P=>D, X=>Q
      real             XP
      data INTEGS / -1 /
c
c++  Code for ARGM is inactive
c      RETURN
c      ENTRY SIVAPE
c++  End
c ********
c START OF CODE
c ********
      IY = 0
      L = NDTF - 1
      do 4680 I = 1, NTE
         INTEG = KORDI
         if (NKDKO .ne. 0) INTEG = KORD(NKDKO + I - 1)
         KQQ = KORD(I + 3)
         K = max(abs(KQQ), 2)
         if (KQQ) 4530, 4520, 4540
 4520    IY = IY + abs(INTEG)
         go to 4670
c ********
c EQUATION IS STIFF, OR IMPLICIT
c ********
 4530    continue
c++  Code for STIFF is inactive
c      KQQ=-KQQ
c      N=KQQ-1
c      JS=abs(KORD(NJSKO+I-1))-1
c      IMPLIC=INTEG
c      INTEG=abs(IMPLIC)-JS
cc.    SET INTEGS FOR STIFF EQUATIONS
c      INTEGS=0
c      IF (K-KSC) 160,160,140
cc.END OF SET-UP FOR STIFF EQUATIONS
c++  End
c ********
c EQUATION IS NOT STIFF
c ********
 4540    N = KQQ
         if (LDT .ne. 0) if (K - KSC) 4570, 4570, 4550
c     DIFFERENCE TABLE HAS NOT BEEN UPDATED
         TP1 = F(I) - F(L + 1)
         if (K - KSC) 4610, 4610, 4590
c END OF SET-UP FOR EQUATIONS WHICH ARE NOT STIFF
c ********
c GET PREDICTED DIFFERENCES FROM UPDATED DIFFERENCE TABLE
c ********
 4550    F(L + K + 1) = F(L + K + 1) * BETA(K + 1)
         TEMP(K) = F(L + K) * BETA(K)
         F(L + K) = TEMP(K)
c LOOP FOR MODIFIED DIVIDED DIFFERENCES
 4560    K = K - 1
         if (K .le. KSC) go to 4580
         TEMP(K) = F(L + K) * BETA(K)
         F(L + K) = TEMP(K) + F(L + K + 1)
         go to 4560
c CODE FOR BACKWARD DIFFERENCES
 4570    F(L + K + 1) = F(L + K + 1)
         TEMP(K) = F(L + K)
         K = K - 1
c
 4580    TEMP(K) = F(L + K)
         F(L + K) = TEMP(K) + F(L + K + 1)
         K = K - 1
         if (K .ne. 0) go to 4580
         go to 4630
c ********
c UPDATE DIFFERENCE TABLE AND GET PREDICTED DIFFERENCES
c ********
c CODE FOR MODIFIED DIVIDED DIFFERENCES
 4590    F(L + K + 1) = (F(L+K+1) + TP1) * BETA(K + 1)
         TEMP(K) = (F(L + K) + TP1) * BETA(K)
         F(L + K) = TEMP(K)
 4600    K = K - 1
         if (K .le. KSC) go to 4620
         TEMP(K) = (F(L + K) + TP1) * BETA(K)
         F(L + K) = TEMP(K) + F(L + K + 1)
         go to 4600
c CODE FOR BACKWARD DIFFERENCES
 4610    F(L + K + 1) = (F(L+K+1) + TP1)
         TEMP(K) = F(L + K) + TP1
         F(L + K) = TEMP(K)
         K = K - 1
c
 4620    TEMP(K) = F(L + K) + TP1
         F(L + K) = TEMP(K) + F(L + K + 1)
         K = K - 1
         if (K .ne. 0) go to 4620
c ********
c COMPUTE Y-S OBTAINED USING INTEGRATION
c ********
c     TEST IF NEXT Y TO BE OBTAINED BY INTERPOLATION
 4630    continue
c++  Code for STIFF is inactive
c      IF (INTEG.EQ.0) GO TO 4662
c++  End
         IY = IY + 1
c     FORM INNER PRODUCT
         XP = C0
         do 4650 J = INTEGS + N + 1, INTEGS + 2, -1
c++  Code for ~{p,x} is active
            XP = XP + G(J, INTEG) * TEMP(J)
c++  Code for {p,x} is inactive
Cc--S Next line special: P=>D, X=>Q
C            XP = XP + real(G(J, INTEG)) * real(TEMP(J))
c++  END
 4650    continue
         K = INTEG + INTEGS
         do 4660 J = K, 1, -1
c++  Code for ~{p,x} is active
            XP = XP + G(1, J) * YN(IY + J)
c++  Code for {p,x} is inactive
Cc--S Next line special: P=>D, X=>Q
C            XP = XP + real(G(1, J)) * real(YN(IY + J))
c++  END
 4660    continue
         Y(IY) = YN(IY) + XP
         INTEG = INTEG - 1
         if (K) 4670, 4670, 4630
c END OF COMPUTING Y-S OBTAINED BY INTEGRATION
c ********
c COMPUTE Y-S OBTAINED USING INTERPOLATION AND DIFFERENTIATION
c ********
c++  Code for STIFF is inactive
cc.    RESTORE INTEGS FOR EQUATIONS WHICH ARE NOT STIFF
c 4662 INTEGS=-1
c      IY=IY+1
cc.    COMPUTE Y USING INTERPOLATION
c      Y(IY)=YN(IY) + F(L+2)
c      IF (KQQ.EQ.1) Y(IY)=YN(IY)
c 4663 INTEG=INTEG+1
c      IF (INTEG.EQ.JS) IF (IMPLIC) 4680,4680,4664
cc.    COMPUTE INTEG-TH DERIVATIVE
c      XP = C0
c 4664 DO 4666 J = KQQ+1, INTEG+1, -1
c         XP = XP + D(J, INTEG) * TEMP(J)
c 4666 CONTINUE
c      IF (INTEG.EQ.JS) GO TO 4667
c      IY=IY+1
c      Y(IY)=XP
c      GO TO 4663
cc.STORE PREDICTED VALUE FOR F
c 4667 CONTINUE
c      F(L+NUMDT)=XP
c++  End
 4670    L = L + NUMDT
 4680    continue
      LDT = -3
      return
      end