C
C
C
C
C
C
C
C
C
C                 *****************************
C                 *****************************
C                 ***                       ***
C                 ***      L A R K I N      ***
C                 ***      (REV. 2.2)       ***
C                 ***                       ***
C                 *****************************
C                 *****************************
C
C
C
C     *********************************************************
C     *                                                       *
C     *            THIS IS REVISION 2.2 OF LARKIN             *
C     *                                                       *
C     *             RELEASE DATE: APRIL 17, '86               *
C     *                                                       *
C     *             ADAPTATION TO IBM-PC: APR. 17, '86        *
C     *                                                       *
C     *********************************************************
C
C
C
C
C
C REFERENCES:
C------------
C
C /1/  P. DEUFLHARD, G. BADER, U. NOWAK:
C      LARKIN - A SOFTWARE PACKAGE FOR THE SIMULATION OF LARGE
C      SYSTEMS ARISING IN CHEMICAL REACTION KINETICS
C      (UNIV. HEIDELBERG, SFB 123: TECHN. REP. 100 (1980))
C      IN:
C      K.H. EBERT, P. DEUFLHARD, W. JAEGER (ED.):
C      MODELLING OF CHEMICAL REACTION SYSTEMS.
C      SPRINGER SERIES CHEM. PHYS. 18 (1981)
C
C /2/  G.BADER, U.NOWAK, P.DEUFLHARD:
C      AN ADVANCED SIMULATION PACKAGE FOR LARGE
C      CHEMICAL REACTION SYSTEMS
C      UNIV. HEIDELBERG, SFB 123: TECHN. REP. 149 (1982)
C
C
C
C
C
C
C---------------------------------------------------------------------
C    DRIVER ROUTINE FOR THE (INTERACTIVE) USE OF THE LARKIN PACKAGE
C---------------------------------------------------------------------
C
C
C
C  BEFORE STARTING THIS DRIVING ROUTINE, PLEASE OBSERVE (AT LEAST)
C  THE FOLLOWING REMARKS:
C
C
C*********************************************************************
C             INTERACTIVE/BATCH USE                                  *
C*********************************************************************
C                                                                    *
C                                                                    *
C  USE LARKIN WITH FULL INTERACTIVE DIALOG FACILITIES:               *
C  SET                                                               *
C    IBATCH = 0  (LINES ENCLOSED IN '**  BATCH **' MAY BE DELETED)   *
C                                                                    *
C  USE LARKIN IN BATCH:                                              *
C  SET                                                               *
C    IBATCH = 1  (LINES ENCLOSED IN '** BATCH **' MAY BE DELETED)    *
C                A DIALOG SESSIN OF LARKIN IS SIMULATED. FILE 'ITERM'*
C                HAS TO BE ALLOCATED TO A DATASET, WHICH CONTAINS    *
C                (LINE BY LINE) THE ANSWERS TO THE QUESTIONS OF THE  *
C                DIALOG SYSTEM. ALL PRINT OUTPUT (INCLUDING DIALOG   *
C                QUESTIONS) IS WRITTEN ON PRINTER.                   *
C                                                                    *
C  OR                                                                *
C    IBATCH = 2  (LINES ENCLOSED IN '**  DIALOG  **' MAY BE DELETED) *
C                 ONLY SUBSYSTEMS 'COMPILE' , 'SIMULATE' AND         *
C                 'NUMERICAL OUTPUT' ARE AVAILABLE. DEFAULT VALUES   *
C                 FOR OUTPUT CONTROL ARE SET INTERNALLY BUT MAY BE   *
C                 CHANGED. ALL PRINT OUTPUT IS WRITTEN ON PRINTER.   *
C                                                                    *
C                                                                    *
C  IMPLEMENTED STANDARD VALUE: IBATCH = 0                            *
C                              ----------                            *
C                                                                    *
C                                                                    *
C*********************************************************************
C             PRINT OUTPUT CONTROL                                   *
C*********************************************************************
C                                                                    *
C  PRINT OUTPUT IS CONTROLLED BY PARAMETER KPRINT. POSSIBLE VALUES   *
C  FOR KPRINT ARE: 0/1/2 (READ 0/1/2 AS: 0 OR 1 OR 2)                *
C  ACCORDING TO KPRINT THE FOLLOWING OUTPUT IS WRITTEN               *
C                                                                    *
C  IN DRIVER:                                                        *
C    KPRINT = 0    ERROR MESSAGES AND DIALOG QUESTIONS               *
C           = 1/2  ERRORS,QUESTIONS AND CPU-TIMES                    *
C                  (ONLY AVAILABLE IF SUITABLE TIME ROUTINES ARE SET *
C                   AND PARAMETER ITIME=1)                           *
C                                                                    *
C  IN SUBSYSTEM 'COMPILE':                                           *
C    KPRINT = 0    ERROR MESSAGES ONLY                               *
C           = 1    ERROR MESSAGES AND WARNINGS                       *
C           = 2    ERROR MESSAGES, WARNINGS AND COMIPLER INPUT       *
C                  LISTING                                           *
C                                                                    *
C  IN SUBSYSTEM 'SIMULATE':                                          *
C    KPRINT = 0    ERROR MESSAGES ONLY                               *
C           = 1    INTEGRATION MONITOR                               *
C           = 2    ENHANCED INTEGRATION MONITOR                      *
C                                                                    *
C  IN SUBSYSTEM 'NUMERICAL OUTPUT':                                  *
C    KPRINT = 0    FIXED STANDARD OUTPUT LISTING                     *
C           = 1,2  DIALOG QUESTIONS AND SELECTED (BY DIALOG) OUTPUT  *
C                                                                    *
C                                                                    *
C  IN CASE OF IBATCH=0 THE OUTPUT IS WRITTEN ON TERMINAL. FOR        *
C  IBATCH=1/2 ON PRINTER.                                            *
C  (I.E.: KPRINT=1 WITH IBATCH=2 IS NOT POSSIBLE FOR SUBSYSTEM 'NUM. *
C   OUTPUT')                                                         *
C                                                                    *
C  FOR IBATCH=0 THE OUTPUT CAN BE DIRECTED TO THE PRINTER IF INSTEAD *
C  OF KPRINT=1/2/3  KPRINT=10/11/12 IS SET.                          *
C                                                                    *
C                                                                    *
C  IMPLEMENTED STANDARD VALUE: KPRINT = 1                            *
C                              ----------                            *
C  (MAY BE TEMPORARY OVERWRITTEN IN DIALOG SESSION)                  *
C                                                                    *
C                                                                    *
C*********************************************************************
C             CPU-TIME ROUTINES                                      *
C*********************************************************************
C                                                                    *
C  THE CALLS OF TIME ROUTINES ARE CONTROLLED BY PARAMETER ITIME      *
C                                                                    *
C  ITIME = 0   NO CALL OF TIME ROUTINES (AND NO ASSOCIATED PRINT     *
C              OUTPUT). TIME ROUTINES MAY BE DUMMY ROUTINES.         *
C                                                                    *
C        = 1   TIME ROUTINES ARE NEEDED TO COMPUTE THE CPU-TIMES     *
C              USED BY EACH SUBSYSTEM (VARIBLE 'STIME') AND OVERALL  *
C              USED CPU-TIME (VARIABLE 'OTIME')                      *
C                                                                    *
C*********************************************************************
C             MAXIMUM DIMENSIONS                                     *
C*********************************************************************
C                                                                    *
C  LARKIN REQUIRES THE DECLARATION OF THE FOLLOWING MAXIMUM          *
C  DIMENSIONS:                                                       *
C  (MAY BE ADAPTED TO MAXIMUM EXPECTED PROBLEM SIZE)                 *
C                                                                    *
C  MAXEL: MAXIMUM NUMBER OF DIFFERENT ELEMENTS USED FOR              *
C         STOICHIOMETRIC CHECK (.LE.10)                              *
C                                                                    *
C  MAXSP: MAXIMUM NUMBER OF DIFFERENT SPECIES IN CHEMICAL REACTION   *
C         SYSTEM (.LT.1000)                                          *
C                                                                    *
C  MAXEQ: MAXIMUM NUMBER OF REACTION EQUATIONS (.LT.1000)            *
C         ( REVERSIBLE EQUATIONS ARE COUNTED BY TWO)                 *
C         ( MAXEQ HAS TO BE SET .GE.MAXSP )                          *
C                                                                    *
C  MAXENT: MAXIMUM NUMBER OF SPECIES ENTRIES IN CHEMICAL EQUATIONS   *
C          (APPROXIMATELY : .LE.4*MAXEQ)                             *
C                                                                    *
C  NRW: DIMENSION OF REAL WORK ARRAY                                 *
C                                                                    *
C                                                                    *
C  NIW: DIMENSION OF INTEGER WORK ARRAY                              *
C                                                                    *
C....................................................................*
C             ACTUAL DIMENSIONS                                      *
C....................................................................*
C                                                                    *
C  THE FOLLOWING DIMENSION ARE ADAPTED TO THE ACTUAL PROBLEM         *
C  AND CALCULATED INTERNALLY:                                        *
C                                                                    *
C                                                                    *
C  IEL: ACTUAL NUMBER OF DIFFERENT ELEMENTS USED FOR                 *
C       BALANCE CHECK                                                *
C                                                                    *
C  NSP: ACTUAL NUMBER OF DIFFERENT SPECIES IN CHEMICAL REACTION      *
C       SYSTEM                                                       *
C                                                                    *
C  NCEQ: ACTUAL NUMBER OF REACTION EQUATIONS                         *
C        ( REVERSIBLE EQUATIONS ARE COUNTED BY TWO)                  *
C                                                                    *
C  LDIM: NUMBER OF NONZEROS IN JACOBIAN                              *
C                                                                    *
C  NLU: NUMBER OF NONZEROS IN DECOMPOSED JACOBIAN                    *
C                                                                    *
C  NLR: NUMBER OF POINTERS REPRESENTING THE STRUCTURE OF CHEMICAL    *
C       SYSTEM AND RIGHT HAND SIDE OF ASSOCIATED ODE-SYSTEM (FCN)    *
C                                                                    *
C*********************************************************************
C             IN/OUTPUT FILES                                        *
C*********************************************************************
C                                                                    *
C  LARKIN REQUIRES THE FOLLOWING FILES:                              *
C  ( FIXED LOGICAL RECORD LENGTH = 80  (CARD IMAGE) )                *
C                                                                    *
C  LOGICAL UNIT               FILE                                   *
C                                                                    *
C      ITIN                   'TERMINAL'                             *
C      ITOUT                  'TERMINAL' (OR 'PRINTER)               *
C      POUT                   'PRINTER'                              *
C                                                                    *
C      CIN                    CHEMICAL INPUT  (DATABASE 'CHEMIN')    *
C                                                                    *
C      DOUT                   NUMERICAL OUTPUT  (DATABASE 'OUT')     *
C                             (FOR STORAGE OF OUTPUT FROM 'SIMULA'   *
C                              USED AS INPUT FOR 'NUM. OUT','PLOT')  *
C                                                                    *
C      HDS                    WORK FILE 1                            *
C                             (FOR INTERMEDIATE OUTPUT STORAGE)      *
C                                                                    *
C      SIN                    SIMULATION INPUT (DATABASE 'INVAL')    *
C                             (FOR STORAGE OF OUTPUT FROM 'CHEMICAL  *
C                              COMPILER', USED AS INPUT FOR 'SIMULA',*
C                              'OPTIONS')                            *
C                                                                    *
C  OPTIONAL:                                                         *
C      INFDAT                 HELP INFORMATON (DATABASE 'INFO')      *
C                                                                    *
C                                                                    *
C  IMPLEMENTED STANDARD VALUES:                                      *
C      ITIN=5 , ITOUT=6 , CIN=8 , HDS=9 , POUT=10 , SIN=11 , DOUT=12 *
C      INFDAT = 20                                                   *
C                                                                    *
C                                                                    *
C*********************************************************************
C             DECLARATIONS                                           *
C*********************************************************************
C                                                                    *
C                                                                    *
C REAL WORK ARRAY:                                                   *
C-----------------                                                   *
C  RW(NRW):  NRW = 12*MAXSP + LDIM + NLU + 2*MAXSP                   *
C                  ( RULE OF THUMB: NRW = (12+7+2*7+2)*MAXSP         *
C                                       = 35*MAXSP                   *
C INTEGER WORK ARRAY:                                                *
C--------------------                                                *
C  IW(NIW):  NIW = MAX( (13*MAXSP + LDIM + NLU + 2*MAXSP) , 360)     *
C                  ( RULE OF THUMB: NIW=36*MAXSP )                   *
C                                                                    *
C                                                                    *
C  COMMON BLOCKS:                                                    *
C----------------                                                    *
C                                                                    *
C  FOR LOGICAL UNITS:                                                *
C     COMMON / LUNIT / ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT              *
C     ( THIS COMMON BLOCK IS USED IN THE WHOLE PACKAGE )             *
C                                                                    *
C                                                                    *
C  FOR ADDITIONAL WORK ARRAYS:                                       *
C                                                                    *
C  GENERAL REMARK:                                                   *
C  ADJUST IN ALL SUBROUTINES THE DIMENSIONING OF ARRAYS DEFINED IN   *
C  COMMON, IF THE COMPILER DOESN'T ACCEPT THE DUMMY DIMENSIONING     *
C  ARRAY(1) INSTEAD OF ARRAY(IDIM) WITH IDIM=EXPLICIT DEFINED        *
C  DIMENSIONS OF ARRAY IN THE DRIVER                                 *
C                                                                    *
C                                                                    *
C     COMMON / LARK1 / LR(MAXENT)                                    *
C                      RULE OF THUMB: MAXENT=4*MAXEQ                 *
C       USED IN SUBROUTINE: SIMULA , FCN , JACOBI                    *
C                                                                    *
C     COMMON / LARK2 / PLR(2*MAXEQ)                                  *
C       USED IN SUBROUTINE: SIMULA , FCN , JACOBI                    *
C                                                                    *
C     COMMON / LARK3 / JA(LDIM)                                      *
C                      RULE OF THUMB: LDIM=7*MAXSP                   *
C       USED IN SUBROUTINE: SIMULA , METASC , JACOBI                 *
C                                                                    *
C     COMMON / LARK4 / IA(MAXSP+1)                                   *
C       USED IN SUBROUTINE: SIMULA , METASC                          *
C                                                                    *
C     COMMON / LARK5 / IDA(MAXSP)                                    *
C       USED IN SUBROUTINE: SIMULA , METASC , JACOBI                 *
C                                                                    *
C     COMMON / KINPAR / RK(MAXEQ)                                    *
C       USED IN SUBROUTINE: SIMULA , SIMOUT , FCN , JACOBI           *
C                                                                    *
C     COMMON / RATE1 / R(L), L=MAX(MAXEQ,90)                         *
C       USED IN SUBROUTINE: SIMULA , METASC , SIMOUT , FCN , JACOBI  *
C                                                                    *
C     COMMON / RATE2 / RINT(MAXEQ)                                   *
C       USED IN SUBROUTINE: SIMULA , METASC , SIMOUT                 *
C                                                                    *
C     COMMON / RATE3 / RTMP(MAXEQ)                                   *
C       USED IN SUBROUTINE: SIMULA , METASC , SIMOUT                 *
C                                                                    *
C                                                                    *
C  FOR MACHINE DEPENDENT CONSTANTS:                                  *
C                                                                    *
C     COMMON / MACHIN / EPMACH , SMALL                               *
C                                                                    *
C      EPMACH: RELATIVE MACHINE PRECISION                            *
C      SMALL: SQUARE ROOT OF SMALLEST POSITIVE MACHINE NUMBER        *
C                                                                    *
C                                                                    *
C  FOR HELP INFORMATION:                                             *
C                                                                    *
C     COMMON / INFORM / IHELP , INFDAT                               *
C                                                                    *
C      IHELP: FLAG                                                   *
C             =0  NO HELP INFORMATION AVAILABLE,                     *
C                 I.E. SUBROUTINE HELP BECOMES A DUMMY ROUTINE AND   *
C                 NO READ FROM LOGICAL UNIT INFDAT OCCURS            *
C             =1  HELP INFORMATION IS WRITTEN IF AN INVALID ANSWER   *
C                 TO A QUESTION OF AN DIALOG BLOCK IS ENTERED        *
C                                                                    *
C      INFDAT: LOGICAL UNIT FOR READ IN SUBROUTINE HELP (FROM THE    *
C              DATABASE CONTAINING THE HELP INFORMATION ('INFO'))    *
C                                                                    *
C                                                                    *
C     COMMON BLOCKS LARK1,..,LARK5,KINPAR,RATE1,RATE2,RATE3 ARE      *
C     USED IN SUBSYSTEM 'SIMULA' ONLY                                *
C     (MAY ALSO BE DECLARED IN SUBROUTINE SIMULA IF OVERLAY TECHNIQUE*
C      IS USED TO REDUCE STORAGE REQUIREMENTS)                       *
C                                                                    *
C                                                                    *
C                                                                    *
C--------------------------------------------------                  *
C  IMPLEMENTED DECLARATIONS HOLDS FOR:                               *
C  MAXEL=10 , MAXSP=400 , MAXEQ=900 , MAXENT=4000 ,                  *
C  RW=14000 , NIW=14400                                              *
C--------------------------------------------------                  *
C                                                                    *
C                                                                    *
C*********************************************************************
C             MACHINE DEPENDENT STATEMENTS                           *
C*********************************************************************
C                                                                    *
C                                                                    *
C  NOTE THAT MACHINE DEPENDENT STATEMENTS ARE ENCLOSED IN ':::'-LINES*
C  THEY APPEAR ONLY IN THIS DRIVER AND IN SUBROUTINE PLOT1.          *
C  CURRENTLY ADAPTED TO IBM 370/168 OF UNIVERSITY OF HEIDELBERG      *
C  AND CALCOMP-SOFTWARE                                              *
C                                                                    *
C                                                                    *
C*********************************************************************
C             ADAPTATION TO IBM-PC                                   *
C*********************************************************************
C                                                                    *
C                                                                    *
C  1.) CHANGE OF TIME ROUTINES IN DRIVER                             *
C  2.) EXPLICIT FILE ALLOCATION IN DRIVER                            *
C  3.) EXPLICIT ENDFILE STATEMENTS IN WHOLE PACKAGE                  *
C  4.) NO CONTINUATION CARDS ALLOWED IN CHEMICAL INPUT               *
C  5.) A-FORMATS FOR WRITE/READ OF VARIABLES DROPPED                 *
C  6.) READ WITH NAMELIST REPLACED                                   *
C  7.) ADDITION OF AN FORTRAN FUNCTION 'DFLOAT'                      *
C  8.) CHANGES IN BLOCK 'OPTION'                                     *
C  9.) CHANGES IN BLOCK 'GRAPH'                                      *
C 10.) CHANGES IN BLOCK 'SIMULA' (ELIMINATE BUG AND CHANGE SCALING)  *
C                                                                    *
C                                                                    *
C*********************************************************************
C
C
C
C  DECLARATIONS
C--------------
C
      DOUBLE PRECISION RW(5000),EPMACH,RK,SMALL,X
CSP       REAL RW(5000),EPMACH,RK,SMALL,X
C
CI4       INTEGER IW(5400)
      INTEGER*2 IW(5400)
C
      REAL FLOAT,R,RINT,RTMP
C
      INTEGER NUM(2)
C
      INTEGER ITERM,ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT
C
CI4       INTEGER TAB(2,10),LR,PLR,JA,IA,IDA
      INTEGER*2 TAB(2,10),LR,PLR,JA,IA,IDA
C
      LOGICAL MASK
C
      COMMON / LUNIT / ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT
C
      COMMON / MACHIN / EPMACH,SMALL
C
      COMMON / INFORM / IHELP,INFDAT
C
      COMMON / LARK1 / LR(4000)
C                      LR(MAXENT)
      COMMON / LARK2 / PLR(1800)
C                      LR(2*MAXEQ)
      COMMON / LARK3 / JA(2800)
C                      JA(LDIM)
      COMMON / LARK4 / IA(401)
C                      IA(MAXSP+1)
      COMMON / LARK5 / IDA(400)
C                      IDA(MAXSP)
      COMMON / KINPAR / RK(900)
C                       RK(MAXEQ)
      COMMON / RATE1 / R(900)
C                      R(L),L:=MAX(MAXEQ,90)
      COMMON / RATE2 / RINT(900)
C                      RINT(MAXEQ)
      COMMON / RATE3 / RTMP(900)
C                      RTMP(MAXEQ)

      COMMON / CCC / ICOUNT
      COMMON / COUNT / NSTEP,NFCN,NANFA,NTFAC,NFAC,NSOL
      COMMON / LOG0 / CLOG0
      COMMON / MA30LE / LP,ABORT1,ABORT2,ABORT3
      COMMON / MA30LF / IRNCP,ICNCP,IRANK,MINIRN,MINICN
      COMMON / MA30LG / DEPS,RMIN
      COMMON / MA30LH / RESID
      COMMON / PLOPC / SPNAME(3,8),XAMI,XAMA,YAMI,YAMA,IIFLAG(10)

c  Description of Integer Common Workspace:
c
c  Name: Common / cigks / igks
c
c  Position      meaning                     value 
c
c
c   1  workstation type (for screen)         set elsewhere 
c   2  connection identifier                 icon=1         
c   3  workstation identifier                iwk=1         
c
c   4  workstation type (for first other wk) set elsewhere
c   5  connection identifier                 icon=2         
c   6  workstation identifier                iwk=2         
c
c   7  workstation type (for first other wk) set elsewhere
c   8  connection identifier                 icon=3         
c   9  workstation identifier                iwk=3         
c
c                                          value     
c
c  10  color index for background            0      
c  11                  forground (white)     1     
c  12                  red                   2    
c  13                  yellow                3   
c  14                  green                 4  
c  15                  cyan                  5 
c  16                  blue                  6      
c  17                  magenta               7     
c  18
c  19
c
c                                          value 
c  20
c  21  aspect source flag no. 1              1
c  22                                        1
c  .
c  .
c  .
c  33  aspect source flag no. 13             1
c
c  49  logical unit for GKS error messages   lgerr=49   
c 
c
      common / cigks / igks(100)
      common / crgks / rgks(100)
c
C
C
C  DEFINE HOLLERITH CONSTANTS (FOR INTERACTIVE DIALOG)
      DATA TAB(1,1)/1HE/ , TAB(2,1)/1He/ ,
     @     TAB(1,2)/1HC/ , TAB(2,2)/1Hc/ ,
     @     TAB(1,3)/1HS/ , TAB(2,3)/1Hs/ ,
     @     TAB(1,4)/1HN/ , TAB(2,4)/1Hn/ ,
     @     TAB(1,5)/1HP/ , TAB(2,5)/1Hp/ ,
     @     TAB(1,6)/1HO/ , TAB(2,6)/1Ho/ ,
     @     TAB(1,7)/1HQ/ , TAB(2,7)/1Hq/
C
C
C
C  DEFINE MAXIMUM DIMENSIONS FOR PROBLEM SIZE
      DATA MAXEL/10/ , MAXSP/200/ , MAXEQ/500/ , MAXENT/4000/
C
C  DEFINE DIMENSIONS OF WORK ARRAYS
      DATA NRW/5000/ , NIW/5400/
C
C
C  OPEN GRAGHIC WINDOW
C---------------------
C
      call gksop0
      call glark0 ('glarkm.init')
      call zlogo2 (48)
C
      call gksdc0(1)
C
C
C  DEFINE LOGICAL UNITS
C----------------------
C
      OPEN ( UNIT=5)
      OPEN ( UNIT=6)
      OPEN ( UNIT=8 , FILE='CHEMIN')
      OPEN ( UNIT=9 , FILE='HDS')
      OPEN ( UNIT=10 , FILE='PRINTER')
      OPEN ( UNIT=11 , FILE='INVAL')
      OPEN ( UNIT=12 , FILE='OUT')
      OPEN ( UNIT=20 , FILE='LINFO')
      ITERM=5
      ITOUTH=6
      CIN=8
      HDS=9
      POUT=10
      SIN=11
      DOUT=12
      INFDAT=20
C
      ITIN=ITERM
      ITOUT=ITOUTH
C
C
C  SET STANDARD VALUE FOR DIALOG/BATCH OPTION OF DRIVER
C------------------------------------------------------
C
      IBATCH=0
C
C
C  SET STANDARD VALUE FOR OUTPUT CONTROL
C---------------------------------------
C
      KPRINT=1
C
C
C  SET STANDARD VALUE FOR CPU-TIME CONTROL
C-----------------------------------------
C
      ITIME=1
C
C
C  SET STANDARD VALUE FOR HELP INFORMATION
C-----------------------------------------
C
      IHELP=1
C
C
C::::::::::::::::::::::::::::::::::::::::::::::::::::
C
C  RELATIVE MACHINE PRECISION
C  (ADAPTED TO IBM 370/168, UNIVERSITY OF HEIDELBERG)
      EPMACH=1.D-16
CSP       EPMACH=1.E-7
C
C  SQUARE-ROOT OF SMALLEST POSITIVE MACHINE-NUMBER
C  (ADAPTED TO IBM 370/168, UNIVERSITY OF HEIDELBERG)
      SMALL=1.D-35
CSP       SMALL=1.E-35
C
      MASK=.TRUE.
C     CALL NOUFL(MASK)
C  (SUPPRESSES STANDARD FIXUP FOR UNDERFLOW)
C::::::::::::::::::::::::::::::::::::::::::::::::::::
C
C  WRITE TOP
C-----------
C
      IF(IBATCH.GE.1) ITOUT=POUT
      WRITE(ITOUT,8000)
8000  FORMAT(///,19X,14H**************,/,19X,14H*   LARKIN   *,/,
     @           19X,14H*  REV. 2.2  *,/,19X,14H**************,/)
C
C
      IF(IBATCH.EQ.2) GOTO 1000
C
C
C  *****  DIALOG  *****  DIALOG  *****  DIALOG *****  DIALOG  *****
C  *****  DIALOG  *****  DIALOG  *****  DIALOG *****  DIALOG  *****
C  *****  DIALOG  *****  DIALOG  *****  DIALOG *****  DIALOG  *****
C
      OTIME=0.
C
C-----------------------------
C  CONTROL INTERACTIVE SESSION
C-----------------------------
C
100   CONTINUE
C
C  RESET STANDARD VALUES
      KFLAG=KPRINT
      ITOUT=ITOUTH
      IF(IBATCH.GE.1) ITOUT=POUT
C
      WRITE(ITOUT,8100)
8100  FORMAT(/,25H INVOKE LARKIN SUBSYSTEM:,/,
     1   54H / EDIT / COMPILE / SIMULATE / NUMERICAL OUTPUT /     ,/,
     2   54H / PLOT / OPTIONS / QUIT / ?                          )
C
C
      CALL ANSNUM (TAB,7,NUM,1,X,IFLAG,1)
C
      IF(NUM(1).GE.0.AND.NUM(1).LE.13) KFLAG=NUM(1)
C
C:::::::::::::::::::::::::::::::::
C     CALL SETTIM
      ITUNIT=100
C  SET COUNTER OF CPU TIME TO ZERO
C:::::::::::::::::::::::::::::::::
C
      GOTO (200,300,400,500,600,700,800) , IFLAG
C
      GOTO 800
C
C  CALL SUBSYSTEM 'EDIT'
C-----------------------
C
200   CONTINUE
C
      CALL EDIT (MAXSP,MAXEQ,MAXEL,KFLAG)
C
C
C::::::::::::::::::::::::::::
C     CALL TIMINT(ITUNIT)
C  TIMINT RETURNS THE NUMBER OF USED CPU-TIME UNITS (UNIT=0.01 SEC)
      STIME=FLOAT(ITUNIT)/100.
C::::::::::::::::::::::::::::
C
      OTIME=OTIME+STIME
      IF(KPRINT.GE.100) WRITE(ITOUT,9100) STIME
102   IF(KFLAG.GE.0) GOTO 100
      WRITE(ITOUT,8200)
8200  FORMAT(//,43H *** ERROR RETURN FROM SUBSYSTEM 'EDIT' ***,/)
      IF(IBATCH.GE.1) GOTO 800
      GOTO 100
C
C  CALL SUBSYSTEM 'CHEMICAL COMPILER'
C--------------------------------------
C
300   CONTINUE
      IF(KFLAG.LT.10) GOTO 301
      KFLAG=KFLAG-10
      ITOUT=POUT
301   CONTINUE
C
      CALL COMPIL (MAXSP,MAXEQ,MAXEL,MAXENT,NIW,NRW,IW,RW,LR,
     @             PLR(1),PLR(MAXEQ+1),KFLAG)
C
      IF(ITIME.EQ.0) GOTO 302
C
C::::::::::::::::::::::::::::
C     CALL TIMINT(ITUNIT)
      STIME=FLOAT(ITUNIT)/100.
C::::::::::::::::::::::::::::
C
      OTIME=OTIME+STIME
      IF(KPRINT.GE.100) WRITE(ITOUT,9100) STIME
302   IF(KFLAG.GE.0) GOTO 100
      WRITE(ITOUT,8300)
8300  FORMAT(//,46H *** ERROR RETURN FROM SUBSYSTEM 'COMPILE' ***,/)
      IF(IBATCH.GE.1) GOTO 800
      GOTO 100
C
C  CALL SUBSYSTEM 'SIMULATION'
C-----------------------------
C
400   CONTINUE
      IF(KFLAG.LT.10) GOTO 401
      KFLAG=KFLAG-10
      ITOUT=POUT
401   CONTINUE
C
      CALL SIMULA (RW,NRW,IW,NIW,KFLAG)
C
      IF(ITIME.EQ.0) GOTO 402
C
C::::::::::::::::::::::::::::
C     CALL TIMINT(ITUNIT)
      STIME=FLOAT(ITUNIT)/100.
C::::::::::::::::::::::::::::
C
      OTIME=OTIME+STIME
      IF(KPRINT.GE.100) WRITE(ITOUT,9100) STIME
402   IF(KFLAG.GE.0) GOTO 100
      WRITE(ITOUT,8400)
8400  FORMAT(//,47H *** ERROR RETURN FROM SUBSYSTEM 'SIMULATE' ***,/)
      IF(IBATCH.GE.1) GOTO 800
      GOTO 100
C
C  CALL SUBSYSTEM 'NUMERICAL OUTPUT'
C-----------------------------------
C
500   CONTINUE
C
      CALL NUMOUT (MAXSP,MAXEQ,MAXEL,NRW,RW,KFLAG)
C
      IF(ITIME.EQ.0) GOTO 502
C
C::::::::::::::::::::::::::::
C     CALL TIMINT(ITUNIT)
      STIME=FLOAT(ITUNIT)/100.
C::::::::::::::::::::::::::::
C
      OTIME=OTIME+STIME
      IF(KPRINT.GE.100) WRITE(ITOUT,9100) STIME
502   IF(KFLAG.GE.0) GOTO 100
      WRITE(ITOUT,8500)
8500  FORMAT(//,43H *** ERROR RETURN FROM SUBSYSTEM 'NUMERICAL,
     @       12H OUTPUT' ***,/)
      IF(IBATCH.GE.1) GOTO 800
      GOTO 100
C
C  CALL SUBSYSTEM 'GRAPHIC OUTPUT'
C---------------------------------
C
600   CONTINUE
      NRW4=NRW*2
CSP       NRW4=NRW
      NIW4=NIW/2
C
      CALL GRAPH (RW,NRW4,IW,NIW4,KFLAG)
C
      IF(ITIME.EQ.0) GOTO 602
C
C::::::::::::::::::::::::::::
C     CALL TIMINT(ITUNIT)
      STIME=FLOAT(ITUNIT)/100.
C::::::::::::::::::::::::::::
C
      OTIME=OTIME+STIME
      IF(KPRINT.GE.100) WRITE(ITOUT,9100) STIME
602   IF(KFLAG.GE.0) GOTO 100
      WRITE(ITOUT,8600)
8600  FORMAT(//,41H *** ERROR RETURN FROM SUBSYSTEM 'GRAPHIC,
     @       12H OUTPUT' ***,/)
      IF(IBATCH.GE.1) GOTO 800
      GOTO 100
C
C  CALL SUBSYSTEM 'OPTIONS'
C--------------------------
C
700   CONTINUE
C
      CALL OPTION (MAXSP,MAXEQ,RW,NRW,IW,NIW,KFLAG)
C
      IF(ITIME.EQ.0) GOTO 702
C
C::::::::::::::::::::::::::::
C     CALL TIMINT(ITUNIT)
      STIME=FLOAT(ITUNIT)/100.
C::::::::::::::::::::::::::::
C
      OTIME=OTIME+STIME
      IF(KPRINT.GE.100) WRITE(ITOUT,9100) STIME
702   IF(KFLAG.GE.0) GOTO 100
      WRITE(ITOUT,8700)
8700  FORMAT(//,46H *** ERROR RETURN FROM SUBSYSTEM 'OPTIONS' ***,/)
      IF(IBATCH.GE.1) GOTO 800
      GOTO 100
C
C  END INTERACTIVE LARKIN SESSION
C--------------------------------
C
800   CONTINUE
      call gkscl0
      WRITE(ITOUT,8800)
8800  FORMAT(///,19X,21HEND OF LARKIN SESSION,///)
      IF(KPRINT.GE.100) WRITE(ITOUT,9200) OTIME
      STOP
C
9100  FORMAT(/,10H CPU-TIME:,F6.2,4H SEC,/)
9200  FORMAT(//,18H OVERALL CPU-TIME:,F6.2,4H SEC,/)
C
C  *****  DIALOG  *****  DIALOG  *****  DIALOG *****  DIALOG  *****
C  *****  DIALOG  *****  DIALOG  *****  DIALOG *****  DIALOG  *****
C
1000  CONTINUE
C
C  *****  BATCH  *****  BATCH  *****  BATCH *****  BATCH  *****
C  *****  BATCH  *****  BATCH  *****  BATCH *****  BATCH  *****
C
C
C  OUTPUT TO PRINTER
      ITOUT=POUT
C
C SET PRINT OPTION FOR COMPILER (KFLAG = 0/1/2)
      KFLAG=2
C
      CALL COMPIL (MAXSP,MAXEQ,MAXEL,MAXENT,NIW,NRW,IW,RW,LR,
     @             PLR(1),PLR(MAXEQ+1),KFLAG)
C
      IF(KFLAG.LT.0) GOTO 1100
C
C SET PRINT OPTION FOR SIMULATION (KFLAG = 0/1/2)
      KFLAG=1
C
      CALL SIMULA (RW,NRW,IW,NIW,KFLAG)
C
      IF(KFLAG.LT.0) GOTO 1100
C
C SET PRINT OPTION FOR NUMERICAL OUTPUT (KFLAG = 0)
      KFLAG=0
C
      CALL NUMOUT (MAXSP,MAXEQ,MAXEL,NRW,RW,KFLAG)
C
      IF(KFLAG.LT.0) GOTO 1100
C
      STOP
C
1100  WRITE(ITOUT,9900)
9900  FORMAT(//,26H LARKIN ENDED DUE TO ERROR,/)
C
      STOP
C  *****  BATCH  *****  BATCH  *****  BATCH *****  BATCH  *****
C  *****  BATCH  *****  BATCH  *****  BATCH *****  BATCH  *****
C  *****  BATCH  *****  BATCH  *****  BATCH *****  BATCH  *****
C
C  END DRIVER ROUTINE FOR LARKIN
C
      END
      DOUBLE PRECISION FUNCTION DFLOAT (I)
C     REAL FUNCTION FLOAT (I)
      DFLOAT=I
C     FLOAT=I
      RETURN
      END
      SUBROUTINE ANSNUM (TAB,NTAB,NUM,NNUM,X,KFLAG,IHELP)
C
C---------------------------------------------------------------------
C  ANSNUM COMPARES AN ANSWER OF THE USER WITH THE FEASIBLE ANSWERS
C  AND SETS AN ASSOCIATED FLAG. OPTIONAL ADDITIONAL GIVEN NUMBERS ARE
C  CHECKED AND CONVERTED FROM STRING IMAGE TO INTEGER (REAL) VALUES
C---------------------------------------------------------------------
C
C
C DATE OF LATEST CHANGE: APRIL 27, '82
C
C
C*********************************************************************
C                                                                    *
C  INPUT PARAMETERS                                                  *
C        TAB:      TABLE OF FEASIBLE ANSWERS                         *
C        NTAB:     NUMBER OF FEASIBLE ANSWERS                        *
C        NNUM:     NUMBER OF (POSSIBLY GIVEN) INTEGER NUMBERS        *
C                  (OR NNUM=-1 IF ONE REAL NUMBER HAS TO BE ANALYSED)*
C        IHELP:    INDICATOR FOR HELP INFORMATION                    *
C                                                                    *
C  OUTPUT PARAMETERS                                                 *
C        NUM(1/2): POSSIBLY ANALYSED INTEGER NUMBERS                 *
C        X:        ANALYSED REAL NUMBER                              *
C                  (ONLY IF NNUM.EQ.-1 AND TABLE SEARCH UNSUCCESSFULL)
C        KFLAG:    ERROR AND TABLE FLAG                              *
C                  .LT.0: INVALID ANSWER                             *
C                  .EQ.0: BLANK LINE                                 *
C                  .GT.0: USER ANSWER FOUND IN TAB(i,KFLAG)          *
C                                                                    *
C*********************************************************************
C
      DOUBLE PRECISION X
CSP       REAL X
CI4       INTEGER TAB(2,NTAB),KEY(80),MKEY,BLANK
      INTEGER*2 TAB(2,NTAB),KEY(80),MKEY,BLANK
      INTEGER NUM(2),IHELP
      INTEGER ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT
C
      COMMON / LUNIT / ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT
C
      DATA BLANK/1H /
C
      K=0
      NUM(1)=-1
      NUM(2)=-1
      DO 1 I=1,73
1     KEY(I)=BLANK
C
2     READ(ITIN,9000) (KEY(L),L=1,73)
      K=K+1
      IF(KEY(73).NE.BLANK) WRITE(ITOUT,8200) (KEY(L),L=63,72)
      KEY(73)=BLANK
      DO 3 I=1,72
      M=I
      IF(KEY(M).NE.BLANK) GOTO 10
3     CONTINUE
      GOTO 20
10    MKEY=KEY(M)
      DO 12 I=1,NTAB
      KFLAG=I
      IF(TAB(1,I).EQ.MKEY) GOTO 30
      IF(TAB(2,I).EQ.MKEY) GOTO 30
12    CONTINUE
      KFLAG=0
      IF(NNUM.NE.-1) GOTO 20
      CALL RNUM (KEY,1,72,X,IFLAG)
      IF(IFLAG.NE.0) KFLAG=-1
      RETURN
C
C  BLANK LINE OR WRONG ANSWER GIVEN
20    CALL HELP (IHELP)
      WRITE(ITOUT,8000) (TAB(1,L),L=1,NTAB)
      IF(K.LT.5) GOTO 2
      WRITE(ITOUT,8100)
      RETURN
C
30    IF(NNUM.EQ.0) RETURN
C
C  ANALYSE NUMBER(S)
40    M=M+1
      IF(M.GT.72) RETURN
      IF(KEY(M).NE.BLANK) GOTO 40
41    M=M+1
      IF(M.GT.72) RETURN
      IF(KEY(M).EQ.BLANK) GOTO 41
      M1=M
      M2=72
42    M=M+1
      IF(M.GT.72) GOTO 43
      IF(KEY(M).NE.BLANK) GOTO 42
      M2=M
43    CALL INUM (KEY,M1,M2,NUM(1),IFLAG)
      IF(IFLAG.NE.0.OR.NNUM.EQ.1) RETURN
      CALL INUM (KEY,M2,72,NUM(2),IFLAG)
      RETURN
C
9000  FORMAT(73A1)
8000  FORMAT(//,24H ENTER NOW YOUR ANSWER: ,/,3H / ,8(A1,3H / ))
8100  FORMAT(//18H TOO MANY ATTEMPTS)
8200  FORMAT(/,27H LINE TRUNCATED AFTER : ...,10A1,/,
     @         24H ENTER CONTINUATION CARD)
C
C  END ANSNUM
C
      END
      SUBROUTINE HELP (I)
C
C---------------------------------------------------------------------
C  HELP LISTS LINES 'IP(I)+1' TO 'IP(I+1)' FROM DATABASE 'INFO', WHICH
C  CONTAINS HELP INFORMATION
C  (FOR IHELP=0 SUBROUTINE HELP IS USED AS A DUMMY ROUTINE; SO THE
C   LOGICAL UNIT 'INFDAT' HAS NOT TO BE ALLOCATED)
C---------------------------------------------------------------------
C
      REAL TXT(18)
      INTEGER ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT
CI4       INTEGER IP(50),CHAR
      INTEGER*2 IP(50),CHAR
C
      COMMON / LUNIT / ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT
C
      COMMON / INFORM / IHELP,INFDAT
C
      IF(IHELP.EQ.0) GOTO 100
C
C
      IF(I.LE.0) GOTO 100
      IF(I.GE.2.AND.I.LT.10) GOTO 100
      IF(I.GT.23) GOTO 100
      IP(1)=1
      IP(2)=61
      IP(10)=61
      IP(11)=75
      IP(12)=108
      IP(13)=145
      IP(14)=152
      IP(15)=156
      IP(16)=161
      IP(17)=167
      IP(18)=171
      IP(19)=176
      IP(20)=197
      IP(21)=201
      IP(22)=206
      IP(23)=210
      IP(24)=274
      IP1=IP(I)
      DO 10 L=1,IP1
10    READ(INFDAT,9000) CHAR
      IP1=IP1+1
      IP2=IP(I+1)
      DO 20 L=IP1,IP2
      READ(INFDAT,9001) (TXT(J),J=1,18)
20    WRITE(ITOUT,9001) (TXT(J),J=1,18)
      REWIND INFDAT
      RETURN
C
100   WRITE(ITOUT,9100)
      RETURN
C
9000  FORMAT(A2)
9001  FORMAT(18A4)
9100  FORMAT(/,3(30H NO HELP INFORMATION AVAILABLE,/))
C
C  END HELP
C
      END
      SUBROUTINE INUM (IN,L1,L2,INTVAL,KFLAG)
C
C---------------------------------------------------------------------
C SUBROUTINE INUM CREATES AN INTEGER VALUE FROM A GIVEN STRING OF DIGITS
C---------------------------------------------------------------------
C
C DATE OF LATEST CHANGE: APRIL 1, '83
C
C*********************************************************************
C                                                                    *
C  INPUT PARAMETERS:                                                 *
C    IN:       INTEGER ARRAY, CONTAINS THE STRING TO BE ANALYSED     *
C                                                                    *
C    L1/L2:    POINTER TO THE FIRST/LAST ELEMENT OF ARRAY IN, WHICH  *
C              HAS TO BE ANALYSED                                    *
C                                                                    *
C                                                                    *
C  OUTPUT PARAMETERS:                                                *
C    INTVAL:   COMPUTED INTEGER VALUE                                *
C                                                                    *
C    KFLAG:    .EQ. 0: SUCCESSFULL ANALYSE                           *
C                   1: BLANK LINE FOUND (NUM:=0)                     *
C                  -1: ANALYSE FAILED                                *
C                                                                    *
C*********************************************************************
C
CI4       INTEGER IN(80),Z(10),BLANK,PLUS,MINUS,CHAR
      INTEGER*2 IN(80),Z(10),BLANK,PLUS,MINUS,CHAR
      INTEGER ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT
      LOGICAL BMINUS
C
      COMMON / LUNIT / ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT
C
      DATA Z(1)/1H0/ ,Z(2)/1H1/ , Z(3)/1H2/, Z(4)/1H3/ , Z(5)/1H4/ ,
     @     Z(6)/1H5/ ,Z(7)/1H6/ , Z(8)/1H7/, Z(9)/1H8/ , Z(10)/1H9/
      DATA BLANK/1H / , PLUS/1H+/ , MINUS/1H-/
C
      BMINUS=.FALSE.
      KFLAG=1
      INTVAL=0
      DO 10 L=L1,L2
      CHAR=IN(L)
      IF(KFLAG.EQ.0) GOTO 1
      IF(CHAR.EQ.BLANK) GOTO 10
      IF(CHAR.EQ.PLUS) GOTO 10
      IF(CHAR.NE.MINUS) GOTO 1
      BMINUS=.TRUE.
      GOTO 10
1     DO 4 K=1,10
      IF(CHAR.NE.Z(K)) GOTO 4
      INTVAL=INTVAL*10 + K-1
      KFLAG=0
      GOTO 10
4     CONTINUE
      IF(BMINUS) INTVAL=-INTVAL
      IF(CHAR.EQ.BLANK) RETURN
      GOTO 100
10    CONTINUE
      IF(BMINUS) INTVAL=-INTVAL
      RETURN
C
C  ERROR HANDLING
100   WRITE(ITOUT,9100) (IN(L),L=L1,L2)
      KFLAG=-1
      RETURN
C
9100  FORMAT(/,14H *** ERROR ***,/,26H INVALID INTEGER DETECTED:,/,
     @       1X,72A1,//)
C
C  END INUM
C
      END
      SUBROUTINE RNUM (IN,L1,L2,X,KFLAG)
C
C---------------------------------------------------------------------
C  SUBROUTINE RNUM CREATES A REAL VALUE FROM A GIVEN STRING OF DIGITS
C---------------------------------------------------------------------
C
C DATE OF LATEST CHANGE: APRIL 27, '82
C
C*********************************************************************
C                                                                    *
C  INPUT PARAMETERS:                                                 *
C  IN :  INTEGER ARRAY, CONTAINS THE STRING TO BE ANALYSED           *
C                                                                    *
C  L1/L2:  POINTER TO THE FIRST/LAST ELEMENT OF ARRAY IN, WHICH      *
C          HAS TO BE ANALYSED                                        *
C                                                                    *
C                                                                    *
C  OUTPUT PARAMETERS:                                                *
C    X   :  COMPUTED REAL VALUE                                      *
C                                                                    *
C    KFLAG:  .EQ. 0: SUCCESSFULL ANALYSE                             *
C                 1: BLANK LINE FOUND (X:=0.)                        *
C                -1: ANALYSE FAILED                                  *
C                                                                    *
C*********************************************************************
C
      DOUBLE PRECISION X,TEN,ZERO,XBASE,DFLOAT
CSP       REAL X,TEN,ZERO,XBASE,FLOAT
      INTEGER ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT
CI4       INTEGER IN(80),Z(10),BLANK,PLUS,MINUS,POINT,CHAR
      INTEGER*2 IN(80),Z(10),BLANK,PLUS,MINUS,POINT,CHAR
CI4       INTEGER DD,EE
      INTEGER*2 DD,EE
CI4       INTEGER DL,EL
      INTEGER*2 DL,EL
      LOGICAL BMINUS,BPOINT,BEXP,BSIGN,BBASE
C
      COMMON / LUNIT / ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT
C
      DATA Z(1)/1H0/ ,Z(2)/1H1/ , Z(3)/1H2/, Z(4)/1H3/ , Z(5)/1H4/ ,
     @     Z(6)/1H5/ ,Z(7)/1H6/ , Z(8)/1H7/, Z(9)/1H8/ , Z(10)/1H9/
      DATA BLANK/1H / , PLUS/1H+/ , MINUS/1H-/ , POINT/1H./
      DATA DD/1HD/ , EE/1HE/ , DL/1Hd/ , EL/1He/
      DATA ZERO/0.D0/ , TEN /1.D1/
CSP       DATA ZERO/0./ , TEN /1.E1/
C
      BMINUS=.FALSE.
      BPOINT=.FALSE.
      BEXP=.FALSE.
      BSIGN=.FALSE.
      BBASE=.FALSE.
      IFLAG=0
      IH=0
      IZ=0
      XBASE=ZERO
      IEXP=0
      X=ZERO
      DO 100 L=L1,L2
      CHAR=IN(L)
      IF(CHAR.EQ.BLANK) GOTO 100
      GOTO (1,5,7) , IFLAG
      IF(CHAR.EQ.POINT) GOTO 10
      IF(CHAR.EQ.PLUS) GOTO 100
      IF(CHAR.EQ.MINUS) GOTO 40
1     IF(CHAR.EQ.DD.OR.CHAR.EQ.EE.OR.CHAR.EQ.DL.OR.CHAR.EQ.EL) GOTO 20
      DO 2 K=1,10
      IF(CHAR.NE.Z(K)) GOTO 2
      XBASE=XBASE*TEN + DFLOAT(K-1)
CSP       XBASE=XBASE*TEN + FLOAT(K-1)
      IF(BPOINT) IEXP=IEXP-1
      BBASE=.TRUE.
      GOTO 100
2     CONTINUE
      GOTO 200
5     IF(CHAR.EQ.MINUS) GOTO 30
      IF(CHAR.EQ.PLUS) GOTO 32
7     DO 9 K=1,10
      IF(CHAR.NE.Z(K)) GOTO 9
      IH=IH*10 + K-1
      IZ=IZ+1
      IF(IZ.LT.3) GOTO 100
      GOTO 300
9     CONTINUE
      GOTO 200
10    BPOINT=.TRUE.
      IFLAG=1
      GOTO 100
20    BEXP=.TRUE.
      IFLAG=2
      GOTO 100
30    BSIGN=.TRUE.
32    IFLAG=3
      GOTO 100
40    BMINUS=.TRUE.
100   CONTINUE
      IF(.NOT.BBASE) GOTO 150
      IF(BEXP.AND.IZ.EQ.0) GOTO 400
      IF(BSIGN) IH=-IH
      X=XBASE*TEN**(IEXP+IH)
      IF(BMINUS) X=-X
      KFLAG=0
      RETURN
C
150   IF(IFLAG.NE.0) GOTO 500
C  BLANK LINE FOUND
98    KFLAG=1
      RETURN
C
C  ERROR HANDLING
200   WRITE(ITOUT,9200) (IN(L),L=L1,L2)
      GOTO 600
300   WRITE(ITOUT,9300) (IN(L),L=L1,L2)
      GOTO 600
400   WRITE(ITOUT,9400) (IN(L),L=L1,L2)
      GOTO 600
500   WRITE(ITOUT,9500) (IN(L),L=L1,L2)
600   KFLAG=-1
      RETURN
C
9200  FORMAT(/14H *** ERROR ***,/,28H NON-ARITHMETIC DIGIT FOUND:,/,
     @  1X,72A1,//)
9300  FORMAT(/14H *** ERROR ***,/,27H TOO MANY DIGITS FOR POWER:,/,
     @  1X,72A1,//)
9400  FORMAT(/14H *** ERROR ***,/,26H MISSING DIGITS FOR POWER:,/,
     @  1X,72A1,//)
9500  FORMAT(/14H *** ERROR ***,/,18H MISSING MANTISSA:,/,1X,72A1,//)
C
C  END RNUM
C
      END
      SUBROUTINE NLCONC (NF,T,CONC,NSP,KFOUND)
C
C*********************************************************************
C
C  SUBROUTINE NLCONC SEARCHES FOR STRING &CON IN FILE 'NF' FROM
C  ACTUAL POSITION TO END OF FILE. IF STRING WAS FOUND,
C  KFOUND IS SET TO +1 AND VALUES OF T AND CONC ARE READ. ELSE
C  KFOUND IS SET TO -1.
C
C*********************************************************************
C
      REAL T,CONC(NSP)
      REAL STRING,CHAR
C
      DATA STRING/4H&CON/
C
10    CONTINUE
      READ(NF,9000,END=20) CHAR
      IF(CHAR.NE.STRING) GOTO 10
C
C  STRING WAS FOUND, READ NOW DATA
      READ(NF,9100) T
      READ(NF,9200) (CONC(L),L=1,NSP)
      KFOUND=1
      RETURN
C  STRING NOT FOUND
20    CONTINUE
      KFOUND=-1
      RETURN
9000  FORMAT(1X,A4)
9100  FORMAT(3X,D21.12)
CSP 9100  FORMAT(3X,E21.12)
9200  FORMAT(5(1X,D13.7,1X))
CSP 9200  FORMAT(5(1X,D13.7,1X))
      END
      SUBROUTINE NLRATE (NF,T,RK,NCEQ,KFOUND)
C
C*********************************************************************
C
C  SUBROUTINE NLRATE SEARCHES FOR STRING &RAT IN FILE 'NF' FROM
C  ACTUAL POSITION TO END OF FILE. IF STRING WAS FOUND,
C  KFOUND IS SET TO +1 AND VALUES OF T AND RK ARE READ. ELSE
C  KFOUND IS SET TO -1.
C
C*********************************************************************
C
      REAL T,RK(NCEQ)
      REAL STRING,CHAR
C
      DATA STRING/4H&RAT/
C
10    CONTINUE
      READ(NF,9000,END=20) CHAR
      IF(CHAR.NE.STRING) GOTO 10
C
C  STRING WAS FOUND, READ NOW DATA
      READ(NF,9100) T
      READ(NF,9200) (RK(L),L=1,NCEQ)
      KFOUND=1
      RETURN
C  STRING NOT FOUND
20    CONTINUE
      KFOUND=-1
      RETURN
9000  FORMAT(1X,A4)
9100  FORMAT(3X,D21.12)
CSP 9100  FORMAT(3X,E21.12)
9200  FORMAT(5(1X,D13.7,1X))
CSP 9200  FORMAT(5(1X,D13.7,1X))
      END
      SUBROUTINE NLIRAT (NF,T,RK,NCEQ,KFOUND)
C
C*********************************************************************
C
C  SUBROUTINE NLIRAT SEARCHES FOR STRING &IRA IN FILE 'NF' FROM
C  ACTUAL POSITION TO END OF FILE. IF STRING WAS FOUND,
C  KFOUND IS SET TO +1 AND VALUES OF T AND RK ARE READ. ELSE
C  KFOUND IS SET TO -1.
C
C*********************************************************************
C
      REAL T,RK(NCEQ)
      REAL STRING,CHAR
C
      DATA STRING/4H&IRA/
C
10    CONTINUE
      READ(NF,9000,END=20) CHAR
      IF(CHAR.NE.STRING) GOTO 10
C
C  STRING WAS FOUND, READ NOW DATA
      READ(NF,9100) T
      READ(NF,9200) (RK(L),L=1,NCEQ)
      KFOUND=1
      RETURN
C  STRING NOT FOUND
20    CONTINUE
      KFOUND=-1
      RETURN
9000  FORMAT(1X,A4)
9100  FORMAT(3X,D21.12)
CSP 9100  FORMAT(3X,E21.12)
9200  FORMAT(5(1X,D13.7,1X))
CSP 9200  FORMAT(5(1X,D13.7,1X))
      END
 
