      PROGRAM DELOAD                                                            
C                                                                               
C WRITTEN BY LIDA A.M. NEJAD,                                                   
C            DEPARTMENT OF MATHEMATICS,                                         
C            UMIST,                                                             
C            MANCHESTER  M60 1QD,                                               
C            UK.                                                                
C                                                                               
C                                      
C                                                                               
C-----------------------------------------------------------------------        
C                                                                               
C    THE  FOLLOWING  IS  A  LISTING  OF  THE PROGRAM DELOAD WHICH WRITES        
C    CHEMICAL RATE EQUATIONS AND APPROPRIATE CONSERVATION EQUATIONS  FOR        
C    A  GIVEN  REACTION  SET.  IN  THE FOLLOWING EXAMPLE, IT SETS UP THE        
C    RELEVANT EQUATIONS FOR AN INTERSTELLAR CLOUD  MODEL  CONSISTING  OF        
C    2411  REACTIONS  AND  247  SPECIES  OF  WHICH  THE ABUNDANCES OF 13        
C    CHEMICAL ELEMENTS, INCLUDING 1 FOR  ELECTRONS,  ARE  CALCULATED  BY        
C    THE  CONSERVATION  EQUATIONS.  (N.B.  THIS  MODEL   INCLUDES   BOTH        
C    NEGATIVE AND POSITIVE IONS).                                               
C                                                                               
C-----------------------------------------------------------------------        
C                                                                               
C INPUT DATA REQUIREMENT                                                        
C ----------------------                                                        
C                                                                               
C (1) A FILE CONTAINING THE REACTION SET.                                       
C (2) A  LIST  CONSISTING  OF  SPECIES  WHOSE CONSERVATION EQUATIONS ARE        
C     REQUIRED TOGETHER WITH PARTICLES SUCH AS  PHOTONS.  THE  CHARACTER        
C     ARRAY  'SPSET'  IS  USED  TO  HOLD  THIS  LIST. A SEPARATE LIST OF        
C     CHEMICAL ELEMENTS, INCLUDING A '+' AND/OR  A  '-'  FOR  IONS;  THE        
C     CHARACTER ARRAY 'ELSET' IS SPECIFIED TO INCORPORATE THIS LIST.            
C                 =====================================                         
C INPUT PROGRAM FLAGS                                                           
C -------------------                                                           
C ISORT   SORTING SPECIFICATIONS; A SPECIES FILE  OR  A  VALID  REACTION        
C         FILE  CAN  BE  PROCESSED  BY  SPECIFYING  ONE OF THE FOLLOWING        
C         VALUES:                                                               
C    =0   IN ADDITION TO (1) AND (2) ABOVE, A DATA FILE  CONTAINING  THE        
C         LIST OF 'N' TIME-DEPENDENT SPECIES IS ALSO GIVEN.  NO  SORTING        
C         IS PERFORMED.                                                         
C    =1   (1),(2)  ARE  PROVIDED, THE  PROGRAM  SELECTS  THE  EXPLICITLY        
C         TIME-DEPENDENT VARIABLE SPECIES AND WRITES THEM TO A FILE.            
C    =2   THE REACTION SET MAY INVOLVE REDUNDANT SPECIES  AND  REACTIONS        
C         OR  IT  HAS  BEEN  MODIFIED  AND  REQUIRES  RENUMBERING.   THE        
C         RELEVENT  REACTIONS  CORRESPONDING TO PRESPECIFIED SPECIES ARE        
C         SELECTED AND A REDUCED, NUMBERED REACTION SET IS CREATED.             
C                                                                               
C IMASS   THE MOLECULAR MASS CALCULATION SPECIFICATION:                         
C    =1   CALCULATION IS REQUIRED                                               
C    =0   CALCULATION IS NOT REQUIRED                                           
C                                                                               
C ISUBS                                                                         
C    =1   THE CHEMISTRY INCLUDES TWO DIGIT SUBSCRIPTS                           
C    =0   THE CHEMISTRY DOES NOT INCLUDE TWO DIGIT SUBSCRIPTS                   
C                                                                               
C INEGE                                                                         
C    =1   NEGATIVE IONS ARE INCLUDED IN THE CHEMISTRY                           
C    =0   NEGATIVE IONS ARE NOT INCLUDED IN THE CHEMISTRY                       
C                 =====================================                         
C INPUT PARAMETERS                                                              
C ----------------                                                              
C NSPEC   NUMBER OF VARIABLES / TIME-DEPENDENT EQUATIONS (WHEN THE EXACT        
C         NUMBER OF SPECIES IS KNOWN)                                           
C MAXSP   > NSPEC IF THE NUMBER OF SPECIES IS NOT KNOWN (I.E.  WHEN  THE        
C         EXACT NUMBER OF SPECIES IS NOT KNOWN)                                 
C NREAC   NUMBER OF REACTIONS                                                   
C MAXRE   >  NREAC  IF  THE NUMBER OF REACTIONS IS NOT KNOWN (I.E.  WHEN        
C         ISORT=2 IS CHOSEN)                                                    
C AMASS   REAL ARRAY OF 'NCONS' ELEMENTS, EACH ELEMENT CONTAINING THE           
C         ATOMIC MASS OF THE CORRESPONDING ELEMENTS IN 'ELSET'                  
C NCONS   NUMBER OF SPECIES CONSERVED / CONSERVATION EQUATIONS                  
C NELEM   NUMBER OF CHEMICAL ELEMENTS AND  COSMIC  PARTICLES  (GIVEN  IN        
C         'SPSET')                                                              
C NONE    NUMBER OF ONE CHARACTER LONG CHEMICAL ELEMENTS INCLUDING A '-'        
C         AND/OR '+' TO ACCOUNT FOR THE IONS (E.G.  NONE=8  FROM 'ELSET'        
C         FOR THE FOLLOWING EXAMPLE PROGRAM)                                    
C NCHAR   MAXIMUM LENGTH OF A SPECIES (IN CHARACTERS)                           
C NHIGH   THE HIGHEST TWO DIGIT SUBSCRIPT PRESENT IN THE CHEMISTRY              
C                                                                               
C A,ALFA, REAL ARRAYS USED FOR COPYING THE CORRESPONDING VALUES  TO  THE        
C BETA    REDUCED RENUMBERED REACTION FILE                                      
C                                                                               
C THE FOLLOWING ARRAYS ARE USED FOR WORKING STORAGE                             
C L,LT,LL,LTL,LP,LTP,SUMM                                                       
C                 =====================================                         
C CHARACTER TYPE VARIABLES                                                      
C ------------------------                                                      
C SPSET  ARRAY  OF  'MAXSP+NELEM'  ELEMENTS  EACH  OF 'NCHAR' CHARACTERS        
C        HOLDING THE SYMBOLIC NAMES OF THE SPECIES                              
C ELSET  ARRAY OF 'NCONS' ELEMENTS EACH  OF  2  CHARACTERS  HOLDING  THE        
C        SYMBOLIC NAMES OF THE CHEMICAL ELEMENTS                                
C RE     ARRAY HOLDING THE SYMBOLIC NAMES OF THE REACTANTS                      
C PR     ARRAY HOLDING THE SYMBOLIC NAMES OF THE PRODUCTS                       
C INDXR  ARRAY; INDEX OF REACTIONS                                              
C INDXS  ARRAY; INDEX OF SPECIES                                                
C INDXC  ARRAY; INDEX OF CONSERVED SPECIES                                      
C SPTMP  ARRAY; TEMPORARY  STORAGE USED WHEN ISORT=2, OF SIZE=5 ALLOWING        
C        FOR TWO REACTANTS AND THREE PRODUCTS                                   
C X      WORKING STORAGE HOLDING TERMS IN A CONSERVATION EQUATION               
C X1     WORKING STORAGE HOLDING NEGATIVE ION TERMS                             
C LOSS   WORKING STORAGE HOLDING THE LOSS TERMS                                 
C PROD   WORKING STORAGE INITIALLY HOLDING THE  PRODUCTION  TERMS,  THEN        
C        THE CONCATENATED PRODUCTION AND LOSS TERMS                             
C SPT    WORKING STORAGE HOLDING A SPECIES FROM THE ARRAY 'SPSET'               
C        THE FOLLOWING VARIABLES REPRESENT A REACTANT:- RE1,RE2                 
C        THE FOLLOWING VARIABLES REPRESENT A PRODUCT:-  PR1,PR2,PR3             
C XR,IXR WORKING STORAGE HOLDING AN ELEMENT FROM THE ARRAY 'INDXR'              
C XS,XN, WORKING STORAGE HOLDING AN ELEMENT FROM THE ARRAY 'INDXS'              
C XJ                                                                            
C XC,IXC WORKING STORAGE HOLDING AN ELEMENT FROM THE ARRAY 'INDXC'              
C                                                                               
C        THE FOLLOWING CHARACTER VARIABLES  ARE  USED AS WORKING STORAGE        
C        AND THEIR SIZES DO NOT NEED TO BE CHANGED:-                            
C        RENUM,MC,MG,MS,NEX1,NEX2,MULT1,MULT2,NPERL                             
C                                                                               
C THE FOLLOWING CHARACTER VARIABLES REPRESENT THE VARIOUS TYPES OF TERMS        
C APPEARING IN THE CONSERVATION, LOSS AND PRODUCTION SUMMATIONS:                
C                                                                               
C TEMP0 = '+Y(I)'                                                               
C TEMP1 = '+M*Y(I)'          RELEVANT TERMS INVOLVED IN CONSERVATION            
C TEMP2 = '+M*Y(I)'          EQUATIONS                                          
C                                                                               
C LTEM1 = '-K(R)*Y(I)'                                                          
C LTEM2 = '-K(R)*Y(I)*D'     VARIOUS FORMS OF POSSIBLE DESTRUCTION TERMS        
C LTEM3 = '-K(R)*Y(I)*X(L)*D'                                                   
C LTEM4 = '-K(R)*Y(I)*Y(J)*D'                                                   
C                                                                               
C PTEM1 = '+K(R)*Y(I)'                                                          
C PTEM2 = '+K(R)*Y(I)*D'                                                        
C PTEM3 = '+K(R)*Y(I)*X(L)*D'                                                   
C PTEM4 = '+K(R)*Y(I)*Y(J)*D'  VARIOUS FORMS OF POSSIBLE FORMATION TERMS        
C PTEM5 = '+K(R)*X(L)'                                                          
C PTEM6 = '+K(R)*X(L)*D'                                                        
C PTEM7 = '+K(R)*X(L)*X(L)*D'                                                   
C                                                                               
C WHERE I=J= 3-DIGIT SPECIES INDEX, L=2-DIGIT CONSERVED  SPECIES  INDEX,        
C R=4-DIGIT  REACTIONS NUMBER INDEX, AND M=1 OR 2-DIGIT FOR COUNTING THE        
C CHEMICAL ELEMENT'S SUBSCRIPT IN A MOLECULE.                                   
C E.G. PTEM3*21, TEMP1*8 AND TEMP2*9                                            
C THE LENGTH OF THE ABOVE PARAMETERS CAN BE MODIFIED IN  CHARACTER  TYPE        
C STATEMENTS.                                                                   
C                                                                               
C N.B.  THE  PRESENCE OF '*D' IN CERTAIN PRODUCTION AND LOSS TERMS ABOVE        
C IS TO ALLOW FOR THE CALCULATION OF FRACTIONAL ABUNDANCES IN MODELS  IN        
C WHICH  DENSITY  VARIES  (E.G.  IN  A  CIRCUMSTELLAR  ENVELOPE   OR   A        
C POST-SHOCK  REGION). IN SUCH CASES THE VALUE OF THE PARAMETER 'D' MUST        
C BE SET TO THE 'DENSITY' AT ANY TIME. HOWEVER, THE VALUE  OF  'D'  MUST        
C BE  SET  EQUAL  TO  1.0  WHEN THE DENSITY IS ASSUMED TO BE CONSTANT IN        
C CERTAIN MODELS. (NOTE THAT THE VALUE OF 'D' (I.E. THE  DENSITY)  NEEDS        
C TO  BE  ASSIGNED  IN  THE  SUBROUTINE DIFFUN, WHEN THE GEAR PACKAGE IS        
C USED; FOR FURTHER  DETAILS  SEE  REFERENCES BELOW).                           
C                                                                               
C-----------------------------------------------------------------------        
C                                                                               
C GENERAL NOTES AND COMMENTS                                                    
C --------------------------                                                    
C                                                                               
C ALTHOUGH  THE PROGRAM ALLOWS FOR UP TO 9999 TWO BODY REACTIONS AND 999        
C SPECIES OF UP TO 7 CHARACTERS EACH, ALL THE ASSOCIATED  VARIABLES  CAN        
C BE  EASILY  CHANGED  TO  THE  DESIRED  VALUES. IT SHOULD BE NOTED THAT        
C THERE ARE LIMITATIONS IN TERMS OF THE FORTRAN COMPILER AVAILABLE.  FOR        
C EXAMPLE,  STANDARD  FORTRAN  DOES  NOT  PERMIT  A  CHARACTER  VARIABLE        
C EXCEEDING 32767 CHARACTERS.  HOWEVER,  THIS  LIMIT  IS  SUFFICIENT  TO        
C SATISFY  THE  REQUIREMENTS FOR A SET OF APPROXIMATELY 10,000 REACTIONS        
C IN WHICH A SPECIES COULD HAVE UP TO ABOUT  2000  PRODUCTION  AND  LOSS        
C TERMS.  ALSO  NOTE  THAT  FOR  LARGE  REACTION  SETS, THE EQUATIONS OF        
C CERTAIN SPECIES ( LIKE  H  AND  C+,  FOR  EXAMPLE)  MAY  CONTAIN  MORE        
C THAN  19  CONTINUATION LINES,  THE  UPPER  LIMIT  ALLOWED BY  STANDARD        
C FORTRAN. IN SUCH  CIRCUMSTANCES,  THE  APPROPRIATE  EQUATIONS  NEED TO        
C BE SPLIT INTO  SUB-EQUATIONS.  DELOAD COPES WITH THE SPLITTING  OF THE        
C VARIABLE SPECIES.  EACH  SPLIT  EQUATION  IS  ASSIGNED TO THE VARIABLE        
C 'YD(J)', WHERE  J=1,..NO.  OF SPLITS, SAY  5  FOR EXAMPLE, AND FINALLY        
C YDOT(I)=YD(1)+...YD(5).   HOWEVER,  SUCH  SPLITTING  FOR  CONSERVATION        
C EQUATIONS LONGER THAN 20 LINES IS NOT YET CARRIED OUT  BY  DELOAD  AND        
C SO NEEDS TO BE DONE MANUALLY USING THE SYSTEM EDITOR.                         
C                                                                               
C                 =====================================                         
C                                                                               
C AN ABRIDGED SAMPLE RECTION SET USED BY DELOAD IS AS FOLLOWS:-                 
C                                                                               
C 1    H       CRP     H+      ELECTR         4.60E-01  0.00      0.0           
C 2    HE      CRP     HE+     ELECTR         5.00E-01  0.00      0.0           
C ...                                                                           
C 57   H+      NH3     NH3+    H              5.20E-09  0.00      0.0           
C 58   H+      C2H2    C2H+    H2             4.30E-09  0.00      0.0           
C 59   H+      H2CO    H2CO+   H              2.00E-09  0.00      0.0           
C ...                                                                           
C 285  C+      C2H4    C3H3+   H              1.02E-09  0.00      0.0           
C 286  C+      CH3OH   CH3+    HCO            2.08E-09  0.00      0.0           
C ...                                                                           
C 999  NH2+    H2O     NH3+    OH             1.00E-10  0.00      0.0           
C 1000 NH2+    H2O     H3O+    NH             2.76E-09  0.00      0.0           
C 1001 NH2+    H2O     NH4+    O              1.45E-10  0.00      0.0           
C ...                                                                           
C 1625 H-      HCN     CN-     H2             3.80E-09  0.00      0.0           
C 1626 C-      O2      O-      CO             4.00E-10  0.00      0.0           
C ...                                                                           
C 1718 HEH+    ELECTR  H       HE             2.00E-07 -0.50      0.0           
C 1719 CH+     ELECTR  C       H              3.00E-07 -0.40      0.0           
C ...                                                                           
C 1944 C       PHOTON  C+      ELECTR         2.16E-10  0.00      2.6           
C 1945 NA      PHOTON  NA+     ELECTR         5.60E-12  0.00      2.0           
C ...                                                                           
C 2078 H+      C-      H       C              2.30E-07 -0.50      0.0           
C ...                                                                           
C 2113 H       CH      C       H2             1.73E-11  0.50   2200.0           
C 2114 H       NH      N       H2             1.73E-11  0.50   2400.0           
C ...                                                                           
C 2411 O       CH      HCO+    ELECTR         2.00E-11  0.44      0.0           
C 9999                                                                          
C                                                                               
C                 =====================================                         
C                                                                               
C AN ABRIDGED SAMPLE SPECIES SET                                                
C                                                                               
C Y(1  )=H        Y(2  )=CO       Y(3  )=H+       Y(4  )=CH   ...               
C Y(6  )=OH       Y(7  )=C2       Y(8  )=NO       Y(9  )=O2   ...               
C ...                                                                           
C ...                                                                           
C Y(226)=NA+      Y(227)=HNO      Y(228)=CCN      Y(229)=C5H  ...               
C Y(231)=S2       Y(232)=O2H      Y(233)=H2O2     Y(234)=CH3O ...               
C                                                                               
C                 =====================================                         
C                                                                               
C AN ABRIDGED SAMPLE OUTPUT OF DELOAD IS AS FOLLOWS:-                           
C                                                                               
C      X( 1) = TOTAL( 1)-(Y(2  )+Y(4  )+2*Y(7  )+Y(12 )+Y(17 )+Y                
C     *        (22 )+Y(25 )+Y(26 )+Y(27 )+2*Y(28 )+3*Y(29 )+Y(30 )+Y(           
C ...                                                                           
C     *        30)+Y(234)+Y(235))                                               
C      X( 2) = TOTAL( 2)-(Y(2  )+Y(6  )+Y(8  )+2*Y(9  )+Y(13 )+Y                
C     *        (19 )+Y(21 )+Y(24 )+Y(25 )+2*Y(30 )+Y(31 )+Y(33 )+Y(37           
C ...                                                                           
C     *        (234))                                                           
C      X( 7) = TOTAL( 7)+(Y(3  )+Y(51 )+Y(64 )+Y(70 )+Y(71 )+Y(7                
C     *        3 )+Y(74 )+Y(75 )+Y(76 )+Y(77 )+Y(78 )+Y(79 )+Y(80 )+Y           
C ...                                                                           
C     *        222)+Y(223)+Y(224)+Y(225)+Y(226))-(+Y(162)+Y(163)+Y(16           
C     *        4)+Y(165)+Y(166)+Y(167))                                         
C ...                                                                           
C      X(13) = TOTAL(13)-(Y(16 )+Y(54 )+Y(55 )+Y(77 )+Y(84 )+Y(1                
C     *        20)+Y(168)+Y(193))                                               
C      YD( 1) =                                                                 
C     *        +K(7   )*X(6 )*D+K(9   )*X(6 )*D+K(9   )*X(6 )*D+K(12            
C ...                                                                           
C     *        )*Y(35 )*D+K(59  )*Y(3  )*Y(37 )*D                               
C      YD( 2) = YD( 1)                                                          
C     *        +K(62  )*Y(3  )*Y(38 )*D+K(63  )*Y(3  )*Y(39 )*D+K(64            
C ...                                                                           
C     *        *Y(5  )*Y(64 )*D+K(220 )*Y(6  )*Y(64 )*D                         
C      YD( 3) = YD( 2)                                                          
C     *        +K(225 )*Y(16 )*Y(64 )*D+K(226 )*Y(11 )*Y(64 )*D+K(235           
C ...                                                                           
C ...                                                                           
C      YD(15) = YD(14)                                                          
C     *        +K(2352)*Y(2  )*Y(15 )*D+K(2399)*Y(28 )*Y(36 )*D+K(240           
C ...                                                                           
C     *        D-K(1686)*Y(163)*D-K(1697)*Y(164)*D))                            
C      YDOT(  1) = +YD(15)                                                      
C     *        +(Y(  1)*(                                                       
C     *        -K(1705)*Y(166)*D-K(1711)*Y(165)*D-K(1715)*Y(167)*D-K(           
C ...                                                                           
C     *        )*Y(6  )*D-K(2403)*Y(9  )*D-K(2404)*Y(24 )*D))                   
C ...                                                                           
C      YDOT( 10) =      +K(1772)*Y(175)*X(7 )*D+K(1774)*Y(122)*X                
C     *        (7 )*D+K(1809)*Y(184)*X(7 )*D+(Y(10 )*(-K(26  )*Y(3  )           
C     *        *D-K(101 )*Y(51 )*D-K(224 )*Y(64 )*D-K(2181)*X(3 )*D-K           
C     *        (2221)*X(2 )*D))                                                 
C      YDOT( 11) =      +K(1757)*Y(119)*X(7 )*D+(Y(11 )*(-K(27                  
C     *        )*Y(3  )*D-K(31  )*Y(3  )*D-K(100 )*Y(51 )*D-K(226 )*Y           
C     *        (64 )*D-K(427 )*Y(76 )*D-K(428 )*Y(76 )*D-K(619 )*Y(83           
C     *         )*D-K(906 )*Y(100)*D-K(1124)*Y(109)*D-K(1301)*Y(125)*           
C     *        D-K(1348)*Y(127)*D-K(1970)-K(2216)*X(2 )*D))                     
C ...                                                                           
C ...                                                                           
C ...                                                                           
C      YDOT(234) =      +K(2335)*Y(6  )*Y(34 )*D+K(2373)*Y(9  )*                
C     *        Y(34 )*D+K(2386)*Y(34 )*Y(63 )*D+(Y(234)*(-K(2147)*Y(1           
C     *          )*D-K(2204)*X(3 )*D-K(2268)*X(2 )*D-K(2342)*Y(6  )*D           
C     *        -K(2377)*Y(9  )*D))                                              
C      YDOT(235) =      +K(1835)*Y(195)*X(7 )*D+K(2286)*Y(22 )*X                
C     *        (4 )*D+(Y(235)*(-K(2201)*X(3 )*D-K(2255)*X(2 )*D))               
C                                                                               
C                 =====================================                         
C                                                                               
C REFERENCES                                                                    
C                                                                               
C BENNETT, A.  1988. IN 'RATE COEFFICIENTS IN ASTROCHEMISTRY' EDS. T.J.         
C        MILLAR AND D.A.WILLIAMS, KLUWER ACADEMIC PUBL. (DORDRECHT), P339.      
C NEJAD, L.A.M. AND MILLAR, T.J. 1987, ASTRON. ASTROPHYS.,183, 279.            
C NEJAD, L.A.M. 1986. PH.D. THESIS, UNIV. OF MANCHESTER.                        
C FARQUHAR, P.R.A. AND MILLAR, T.J., 1993, CCP7 NEWSLETTER NO. 18, 6                                                                              
C-----------------------------------------------------------------------        
C                                                                               
      PARAMETER(MAXRE=4000,MAXSP=394,NCONS=2,NELEM=5)                           
      PARAMETER(ISORT=2,IMASS=0,ISUBS=1,INEGE=1)                                
      PARAMETER(NPERL=54,IPL=NPERL-1,LCHAR=18*NPERL,LCHAR1=LCHAR+40,            
     *          INUM1=NELEM-1,NHIGH=12,NTWO=NHIGH-9,NTOT=MAXSP+NELEM)           
      PARAMETER(NONE=2,NCHAR=7)                                                 
      CHARACTER*1  NEX1,MULT1,RENUM(0:9)                                        
      CHARACTER*3  XC1,ELSET(INUM1),SM(INUM1),MS(NTWO),INDXC(INUM1),XC,
     *             NEX2,MULT2                                                   
      CHARACTER*3  XS,XN,XJ,INDXS(MAXSP),PR3,PR4
      CHARACTER*4  XR,IXR,INDXR(MAXRE)                                          
      CHARACTER*7  RE(2,MAXRE),SPECI(MAXSP),PR(4,MAXRE),RE1,RE2,PR1,            
     *             PR2,SPSET(NTOT),SPTMP(6),SPT,DEL*9                   
      CHARACTER    MG*10,MC*9,WS*200000,                                          
     *             TEMP1*9,TEMP0*7,TEMP2*10,                                    
     *             LTEM1*9,LTEM2*18,LTEM3*20,LTEM4*18,                          
     *             PTEM1*16,PTEM2*25,PTEM3*26,PTEM4*15,PTEM5*23,                
     *             X1*2000,X*9000,LOSS*200000,PROD*200000                          
      REAL AMASS(NCONS),SUMM(MAXSP),A(MAXRE),ALFA(MAXRE),BETA(MAXRE)            
      INTEGER L(NCONS),LT(NCONS),LL(MAXSP),LTL(MAXSP),LTP(MAXSP),               
     *        LP(MAXSP)                                                         
      DATA (SPSET(I),I=1,NELEM)/'ELECTR',                                  
     *     'CRP','PHOTON','CRPHOT','     '/,                    
     *     (RENUM(I),I=0,9)/'0','1','2','3','4','5','6','7','8','9'/,           
     *     MC/'23456789'/, MG/'123456789'/,
     *     (ELSET(I),I=1,NCONS)/'-','+'/
                                  
C                                                                               
C OPEN APPROPRIATE DATA AND OUTPUT FILES                                        
C SPECIES FILE IS ON UNIT 10                                                    
C REACTION FILE IS ON UNIT 11                                                   
C OUTPUT FILE IS ON UNIT 12 
C NEW RATE FILE IS ON UNIT 13                                                    
C N.B. THE CONTENTS OF THE OPEN STATEMENTS CAN BE MACHINE DEPENDENT             
      OPEN(UNIT=10,FILE="species94.data")                           
      OPEN(UNIT=11,FILE="rate95.data")
      OPEN(UNIT=12,FILE="odes.f")
      OPEN(UNIT=13,FILE="newrate.data")                                                        
C******************************* SECTION 1 *****************************        
      NSPEC = MAXSP                                                             
      IF (ISORT.EQ.0) THEN                                                      
C ----------------------------- (ISORT=0) ------------------------------        
         READ(10,1002)(INDXS(J),SPECI(J),J=1,NSPEC)                             
         WRITE(12,1003)(INDXS(J),SPECI(J),J=1,NSPEC)                            
         DO 15 J=1,NSPEC                                                        
         SPSET(NELEM+J) = SPECI(J)                                              
 15      CONTINUE                                                               
         READ(11,1001)                                                          
         DO 35 J=1,MAXRE                                                        
         READ(11,1004)INDXR(J),(RE(I,J),I=1,2),(PR(I,J),I=1,4),AJ,ALFAJ,        
     *             BETAJ,DEL                                                          
         IF(INDXR(J).EQ.'9999') GO TO 40                                       
 35      CONTINUE                                                               
 40      NREAC=MAXRE                                                            
      ELSE IF (ISORT.EQ.1) THEN                                                 
C ----------------------------------------------------------------------        
C ----------------------------- (ISORT=1) ------------------------------        
         READ(11,1001)                                                          
         DO 125 J=1,MAXRE                                                       
         READ(11,1004)INDXR(J),(RE(I,J),I=1,2),(PR(I,J),I=1,4),AJ,ALFAJ,        
     *             BETAJ,DEL         
         WRITE(12,1004)INDXR(J),(RE(I,J),I=1,2),(PR(I,J),I=1,4),AJ,             
     *             ALFAJ,BETAJ,DEL                                              
         IF (INDXR(J).EQ.'9999') GO TO 130                                      
 125     CONTINUE                                                               
 130     NREAC = J-1                                                            
         DO 140 I=1,MAXSP                                                       
 140     INDXS(I) = INDXR(I)                                                    
         NCODE=NELEM                                                            
         DO 180 IR=1,NREAC                                                      
         DO 170 M=1,2                                                           
         DO 150 N1=1,NCODE+1                                                    
         IF (RE(M,IR).EQ.SPSET(N1)) GO TO 160                                   
         IF (N1.GT.NCODE) THEN                                                  
            NCODE = NCODE+1                                                     
            NSPEC = NCODE-NELEM                                                 
            SPECI(NSPEC) = RE(M,IR)                                             
            SPSET(NCODE) = RE(M,IR)                                             
            GO TO 160                                                           
         END IF                                                                 
 150     CONTINUE                                                               
 160     GO TO (170,180)M                                                       
 170     CONTINUE                                                               
 180     CONTINUE                                                               
         WRITE(12,1003)(INDXS(J),SPECI(J),J=1,NSPEC)                            
      ELSE IF (ISORT.EQ.2) THEN                                                 
C ----------------------------------------------------------------------        
C ----------------------------- (ISORT=2) ------------------------------        
         READ(10,1002)(INDXS(J),SPECI(J),J=1,NSPEC)                             
         WRITE(12,1003)(INDXS(J),SPECI(J),J=1,NSPEC)                            
         DO 210 J=1,NSPEC                                                       
         SPSET(NELEM+J) = SPECI(J)                                              
 210     CONTINUE                                                               
         IRUSED = 0                                                             
         READ(11,1001)                                                          
         DO 275 J=1,MAXRE                                                       
         READ(11,1004)IXR,RE1,RE2,PR1,PR2,PR3,PR4,AJ,ALFAJ,BETAJ,DEL
         SPTMP(1) = RE1                                                         
         SPTMP(2) = RE2                                                         
         SPTMP(3) = PR1                                                         
         SPTMP(4) = PR2                                                         
         SPTMP(5) = PR3 
         SPTMP(6) = PR4                                                        
         DO 255 M=1,6                                                           
         DO 245 K1=1,NTOT                                                       
         IF (SPTMP(M).EQ.SPSET(K1)) THEN                                        
            GO TO (255,255,255,255,255,260)M                                        
         END IF                                                                 
 245     CONTINUE                                                               
         GO TO 275                                                              
 255     CONTINUE                                                               
 260     IRUSED = IRUSED+1                                                      
         RE(1,IRUSED) = RE1                                                     
         RE(2,IRUSED) = RE2                                                     
         PR(1,IRUSED) = PR1                                                     
         PR(2,IRUSED) = PR2                                                     
         PR(3,IRUSED) = PR3
         PR(4,IRUSED) = PR4                                                     
         A(IRUSED) = AJ                                                         
         ALFA(IRUSED) = ALFAJ                                                   
         BETA(IRUSED) = BETAJ                                                   
                                                                                
C NUMBER VALID REACTIONS                                                        
                                                                                
         IF (IRUSED.LT.10) THEN                                                 
            INDXR(IRUSED) = RENUM(IRUSED)                                       
         ELSE IF (IRUSED.LT.100) THEN                                           
            J1 = INT(IRUSED/10)                                                 
            J2 = IRUSED-J1*10                                                   
            INDXR(IRUSED) = RENUM(J1)//RENUM(J2)                                
         ELSE IF (IRUSED.LT.1000) THEN                                          
            J1 = INT(IRUSED/100)                                                
            IR1 = IRUSED-J1*100                                                 
            J2 = INT(IR1/10)                                                    
            J3 = IR1-J2*10                                                      
            INDXR(IRUSED) = RENUM(J1)//RENUM(J2)//RENUM(J3)                     
         ELSE                                                                   
            J1 = INT(IRUSED/1000)                                               
            IR1 = IRUSED-J1*1000                                                
            J2 = INT(IR1/100)                                                   
            IR2 = IR1-J2*100                                                    
            J3 = INT(IR2/10)                                                    
            IR3 = IR2-J3*10                                                     
            J4 = IR3                                                            
            INDXR(IRUSED)= RENUM(J1)//RENUM(J2)//RENUM(J3)//RENUM(J4)           
         END IF                                                                 
         IF (IXR.EQ.'9999') GO TO 280                                           
 275     CONTINUE                                                               
 280     NREAC = IRUSED                                                         
                                                                                
C  WRITE VALID REACTION FILE TO UNIT 13                                 
                                                                                                                      
         WRITE(13,1001)                                                         
         DO 285 J=1,IRUSED                                                      
         WRITE(13,1004)INDXR(J),(RE(I,J),I=1,2),(PR(I,J),I=1,4),                
     *              A(J),ALFA(J),BETA(J)                                        
 285     CONTINUE                                                               
         CLOSE(UNIT=13)                                                         
      END IF                                                                    
                                                                                
C ----------------------------------------------------------------------        
                                                                                
C LOAD ARRAY MS                                                                 
                                                                                
      DO 290 I=1,NTWO                                                           
      I2 = 9+I                                                                  
      MS(I) = INDXS(I2)                                                         
 290  CONTINUE                                                                  
                                                                                
C ------------------- CALCULATION OF MOLECULAR MASSES ------------------        
C IF THIS CALCULATION IS NOT REQUIRED, THE FOLLOWING LINES UP TO                
C STATEMENT 400 AND,ARRAYS SUMM,AMASS,MG AND FLAG 'IMASS' MAY BE                
C REMOVED FROM THE CODE.                                                        
C ----------------------------------------------------------------------        
      IF (IMASS.EQ.1) THEN                                                      
         DO 305 I=1,INUM1                                                       
         SM(I) = ELSET(I)                                                       
 305     CONTINUE                                                               
         NMASS = NCONS                                                          
         DO 380 NM = 1,NSPEC                                                    
         SPT = SPECI(NM)                                                        
         I = 0                                                                  
         I1 = 0                                                                 
         J1 = 0                                                                 
         J2 = 0                                                                 
  310    I = I+1+J1+J2                                                          
         J1 = 0                                                                 
         J2 = 0                                                                 
         I1 = 0                                                                 
         IF (SPT(I:I).NE.' ') THEN                                              
            DO 360 K=1,NMASS                                                    
            DO 315 M=NONE+1,NCONS                                               
            IF (I.GT.NCHAR-1) GO TO 320                                         
            IF (SPT(I:I+1).EQ.SM(M)) THEN                                       
               I1 = 1                                                           
               J1 = 1                                                           
               GO TO 320                                                        
            END IF                                                              
 315        CONTINUE                                                            
 320        IF ((I+I1).LE.NCHAR) THEN                                           
               IF (K.GT.NONE) I1=1                                              
               IF (SPT(I:I+I1).EQ.SM(K)(1:1+I1)) THEN                           
C                                                                               
C CHECK FOR TWO DIGIT ELEMENTAL SUBSCRIPTS AND TAKE THEM INTO ACCOUNT           
C                                                                               
               IF (ISUBS.NE.0) THEN                                             
                  NEX2 = SPT(I+1+I1:I+2+I1)                                     
                  DO 325 J=1,NTWO                                               
                  J3 = 9+J                                                      
                  IF (NEX2.EQ.MS(J)) THEN                                       
                     MJ = J3                                                    
                     J2 = 2                                                     
                     GO TO 350                                                  
                  END IF                                                        
 325              CONTINUE                                                      
               END IF                                                           
C                                                                               
C CHECK FOR ONE DIGIT ELEMENTAL SUBSCRIPTS                                      
C                                                                               
               NEX1 = SPT(I+1+I1:I+1+I1)                                        
               IF(.NOT.(NEX1.EQ.' '.OR.NEX1.EQ.'+'.OR.NEX1.EQ.'-')) THEN        
                  DO 335 J=1,9                                                  
                  IF (NEX1.EQ.MG(J:J)) THEN                                     
                     MJ =J                                                      
                     J2 = 1                                                     
                     GO TO 350                                                  
                  END IF                                                        
 335              CONTINUE                                                      
               END IF                                                           
               MJ = 1                                                           
 350           SUMM(NM) = SUMM(NM)+AMASS(K)*MJ                                  
               GO TO 310                                                        
               END IF                                                           
            END IF                                                              
 360        CONTINUE                                                            
            GO TO 310                                                           
         END IF                                                                 
 380     CONTINUE                                                               
C                                                                               
         WRITE(12,1001)                                                         
         WRITE(12,390)(SUMM(J),J=1,NSPEC)                                       
 390     FORMAT(5(3X,F6.2,1X))                                                  
      END IF                                                                    
C-----------------------------------------------------------------------        
 400  CONTINUE                                                                  
C****************************** SECTION 2 ******************************        
                                                                                
C------------------- SET UP CONSERVATION EQUATIONS ---------------------        
                                                                                
      IFLAG = 0                                                                 
      LTMC0 = LEN(TEMP0)                                                        
      LTMC1 = LEN(TEMP1)                                                        
      LTMC2 = LEN(TEMP2)                                                        
      DO 500 K=1,NCONS                                                          
      L(K) = LT(K)                                                              
      DO 460 N1 = 1,NSPEC                                                       
      SPT = SPECI(N1)                                                           
      XS = INDXS(N1)                                                            
                                                                                
      I = 0                                                                     
      I1 = 0                                                                    
      J2 = 0                                                                    
      J1 = 0                                                                    
  410 I = I+1+J1+J2                                                             
      I1 = 0                                                                    
      J2 = 0                                                                    
      J1 = 0                                                                    
      IF (SPT(I:I).NE.' ') THEN                                                 
         DO 415 M=NONE+1,NCONS                                                  
         IF (I.GT.NCHAR-1) GO TO 420                                            
         IF (SPT(I:I+1).EQ.ELSET(M)) THEN                                       
            I1 = 1                                                              
            J1 = 1                                                              
            GO TO 420                                                           
         END IF                                                                 
 415     CONTINUE                                                               
 420     IF ((I+I1).LE.NCHAR) THEN                                              
            IF (K.GT.NONE) I1=1                                                 
            IF (SPT(I:I+I1).EQ.ELSET(K)(1:1+I1)) THEN                           
C                                                                               
C CHECK FOR TWO DIGIT ELEMENTAL SUBSCRIPTS AND TAKE THEM INTO ACCOUNT           
C                                                                               
               IF (ISUBS.NE.0) THEN                                             
                  NEX2 = SPT(I+1+I1:I+2+I1)                                     
                  DO 430 J=1,NTWO                                               
                  IF (NEX2.EQ.MS(J)) THEN                                       
                     MULT2 = NEX2                                               
                     TEMP2 = '+'//MULT2//'*Y('//XS//')'                         
                     L2 = LTMC2+L(K)                                            
                     LT(K) = L2                                                 
                     L1 = L(K)+1                                                
                     WS(L1:L2) = TEMP2                                          
                     J2 = 2                                                     
                     GO TO 450                                                  
                  END IF                                                        
 430              CONTINUE                                                      
               END IF                                                           
C                                                                               
C CHECK FOR ONE DIGIT ELEMENTAL SUBSCRIPTS                                      
C                                                                               
               NEX1 = SPT(I+1+I1:I+1+I1)                                        
               IF(.NOT.(NEX1.EQ.' '.OR.NEX1.EQ.'+'.OR.NEX1.EQ.'-')) THEN        
                  DO 440 J=1,9                                                  
                  IF (NEX1.EQ.MC(J:J)) THEN                                     
                     MULT1 = MC(J:J)                                            
                     TEMP1 = '+'//MULT1//'*Y('//XS//')'                         
                     L2 = LTMC1+L(K)                                            
                     LT(K) = L2                                                 
                     L1 = L(K)+1                                                
                     WS(L1:L2) = TEMP1                                          
                     J2 = 1                                                     
                     GO TO 450                                                  
                  END IF                                                        
 440              CONTINUE                                                      
               END IF                                                           
               TEMP0 = '+'//'Y('//XS//')'                                       
               L2 = LTMC0+L(K)                                                  
               LT(K) = L2                                                       
               L1 = L(K) +1                                                     
               WS(L1:L2) = TEMP0                                                
            END IF                                                              
         END IF                                                                 
 450     L(K) = LT(K)                                                           
         GO TO 410                                                              
      END IF                                                                    
 460  CONTINUE                                                                  
      ID = LT(K)+1                                                              
      WS(ID:ID+60) = ' '                                                        
      IF (INEGE.EQ.1) THEN                                                      
         IF (ELSET(K).EQ.'-') THEN                                              
            LCON1 = LT(K)+3                                                     
            X1 = '-('//WS(1:LT(K))//')'                                         
            GO TO 500                                                           
         END IF                                                                 
         IF (ELSET(K).EQ.'+') THEN                                              
            LCON = LCON1+LT(K)+2                                                
            X = '+('//WS(2:LT(K))//')'//X1                                      
            IFLAG = 1                                                           
            IB = K-1                                                            
            IC = IB                                                             
            GO TO 480                                                           
         END IF                                                                 
         IF (IFLAG.EQ.1) THEN                                                   
            IB = K-1                                                            
            IC = IB                                                             
            GO TO 475                                                           
         END IF                                                                 
      END IF                                                                    
      IB = K                                                                    
      IC = IB                                                                   
 475  X = '-('//WS(2:LT(K))//')'                                                
      LCON = LT(K)+2                                                            
 480  IF (LCON.GT.40) THEN                                                      
         WRITE(12,1005)IC,IB,X(:40)                                             
         NLINE = (LCON-40)/NPERL                                                
         IF ((LCON-40).NE.NLINE*NPERL) NLINE=NLINE+1                            
         DO 495 J=1,NLINE                                                       
         J1 = J*40+(4*(J-1)+1)+(J-1)*10                                         
         J2 = J1+IPL                                                            
         WRITE(12,1006)X(J1:J2)                                                 
 495     CONTINUE                                                               
      ELSE                                                                      
         WRITE(12,1005)IC,IB,X(:40)                                             
      END IF                                                                    
 500  CONTINUE                                                                  
C****************************** SECTION 3 ******************************        
                                                                                
      DO 505 J=1,INUM1                                                         
      INDXC(J) = INDXS(J)                                                       
  505 CONTINUE                                                                  
C                                                                               
      LTML1 = LEN(LTEM1)                                                        
      LTML2 = LEN(LTEM2)                                                        
      LTML3 = LEN(LTEM3)                                                        
      LTML4 = LEN(LTEM4)                                                        
C     LTML5 = LEN(LTEM5)                                                        
      LTMP1 = LEN(PTEM1)                                                        
      LTMP2 = LEN(PTEM2)                                                        
      LTMP3 = LEN(PTEM3)                                                        
      LTMP4 = LEN(PTEM4)                                                        
      LTMP5 = LEN(PTEM5)                                                        
C     LTMP6 = LEN(PTEM6)                                                        
C     LTMP7 = LEN(PTEM7)                                                        
      IBCK = 5+LEN(XS)                                                          
C-------------------------- SET UP LOSS TERMS --------------------------        
C                                                                               
      DO 1000 NS=1,NSPEC                                                        
      XS = INDXS(NS)                                                            
      DO 550 NR= 1,NREAC                                                        
      XR = INDXR(NR)                                                            
      LL(NS) = LTL(NS)                                                          
      DO 510 M=1,2                                                              
      IF (SPECI(NS).EQ.RE(M,NR)) THEN                                           
         IF (M.EQ.1) THEN                                                       
            M1=M+1                                                              
         ELSE                                                                   
            M1=M-1                                                              
         END IF                                                                 
         GOTO 530                                                               
      END IF                                                                    
 510  CONTINUE                                                                  
      GO TO 550                                                                 
 530  DO 540 IC=1,INUM1                                                         
      XC = INDXC(IC)                                                            
      IF (RE(M1,NR).EQ.SPSET(IC)) THEN                                          
         IF (SPSET(IC).EQ.'PHOTON'.OR.SPSET(IC).EQ.'CRP'
     *       .OR.SPSET(IC).EQ.'CRPHOT') THEN                  
            LTEM1 = '-K('//XR//')'                                              
            L2 = LTML1+LL(NS)                                                   
            LTL(NS) =L2                                                         
            L1 = LL(NS)+1                                                       
            LOSS(L1:L2) = LTEM1                                                 
C        ELSE IF (SPSET(IC).EQ.'   ') THEN                                      
C           LTEM2 = '-K('//XR//')*D'                                            
C           L2 = LTML2+LL(NS)                                                   
C           LTL(NS) =L2                                                         
C           L1 = LL(NS)+1                                                       
C           LOSS(L1:L2) = LTEM2                                                 
         ELSE                                                                   
            LTEM2 = '-K('//XR//')*X('//XC//')*D'                                
            L2 = LTML2+LL(NS)                                                   
            LTL(NS) = L2                                                        
            L1 = LL(NS)+1                                                       
            LOSS(L1:L2) = LTEM2                                                 
         END IF                                                                 
         GOTO 550                                                               
      END IF                                                                    
 540  CONTINUE                                                                  
      DO 545 J=1,NSPEC                                                          
      XJ = INDXS(J)                                                             
      IF (RE(M1,NR).EQ.RE(M,NR).AND.RE(M1,NR).EQ.SPECI(J)) THEN                 
         LTEM3 = '-2*K('//XR//')*Y('//XJ//')*D'                                 
         L2 = LTML3+LL(NS)                                                      
         LTL(NS) = L2                                                           
         L1 = LL(NS)+1                                                          
         LOSS(L1:L2) = LTEM3                                                    
         GO TO 550                                                              
      ELSE IF (RE(M1,NR).EQ.SPECI(J)) THEN                                      
         LTEM4 = '-K('//XR//')*Y('//XJ//')*D'                                   
         L2 = LTML4+LL(NS)                                                      
         LTL(NS) = L2                                                           
         L1 = LL(NS)+1                                                          
         LOSS(L1:L2) = LTEM4                                                    
         GO TO 550                                                              
      END IF                                                                    
 545  CONTINUE                                                                  
 550  CONTINUE                                                                  
C-----------------------------------------------------------------------        
C WRITE LOSS TERMS TO OUTPUT FILE AS A CHECK ON FAILURE - THIS WILL ONLY        
C BE PERFORMED IF THE 'C' IN COLUMN 1 IS REPLACED WITH A BLANK IN THE           
C LINES UP TO STATEMENT NUMBER 590.                                             
C                                                                               
C      IF (LTL(NS).GT.40) THEN                                                  
C      WRITE(12,570)NS,LOSS(:40)                                                
C 570 FORMAT(6X,'LOSS(',I3,')= ',A40)                                           
C      NLINE = (LTL(NS)-40)/NPERL                                               
C      IF ((LTL(NS)-40).GT.NLINE*NPERL) NLINE = NLINE+1                         
C      DO 580 J=1,NLINE                                                         
C      J1 = J*40+(4*(J-1)+1)+(J-1)*10                                           
C      J2 = J1+NPERL-1                                                          
C      WRITE(12,1006)LOSS(J1:J2)                                                
C 580  CONTINUE                                                                 
C      ELSE                                                                     
C      WRITE(12,570) NS,LOSS(:40)                                               
C 590  END IF                                                                   
C-----------------------------------------------------------------------        
C----------------------- SET UP PRODUCTION TERMS -----------------------        
C                                                                               
      DO 820 NR=1,NREAC                                                         
      XR = INDXR(NR)                                                            
      LP(NS) = LTP(NS)                                                          
      DO 815 MP=1,3                                                             
      IF (SPECI(NS).EQ.PR(MP,NR)) THEN                                          
         DO 750 N1=1,NSPEC                                                      
         XN = INDXS(N1)                                                         
         DO 640 M=1,2                                                           
         IF (RE(M,NR).EQ.SPECI(N1)) THEN                                        
            IF (M.EQ.1) THEN                                                    
               M1=M+1                                                           
            ELSE                                                                
               M1=M-1                                                           
            END IF                                                              
            GOTO 680                                                            
         END IF                                                                 
 640     CONTINUE                                                               
         GOTO 750                                                               
 680     DO 690 IC=1,INUM1                                                      
         XC = INDXC(IC)                                                         
         IF (RE(M1,NR).EQ.SPSET(IC)) THEN                                       
             IF (SPSET(IC).EQ.'PHOTON'.OR.SPSET(IC).EQ.'CRP'
     *        .OR.SPSET(IC).EQ.'CRPHOT') THEN              
               PTEM1 = '+K('//XR//')*Y('//XN//')'                               
               L2 = LTMP1+LTP(NS)                                               
               LTP(NS) = L2                                                     
               L1 = LP(NS)+1                                                    
               PROD(L1:L2) = PTEM1                                              
C           ELSE IF (SPSET(IC).EQ.'  ') THEN                                    
C              PTEM2 = '+K('//XR//')*Y('//XN//')*D'                             
C              L2 = LTMP2+LTP(NS)                                               
C              LTP(NS) = L2                                                     
C              L1 = LP(NS)+1                                                    
C              PROD(L1:L2) = PTEM2                                              
            ELSE                                                                
               PTEM2 = '+K('//XR//')*Y('//XN//')*X('//XC//')*D'                 
               L2 = LTMP2+LTP(NS)                                               
               LTP(NS) = L2                                                     
               L1 = LP(NS)+1                                                    
               PROD(L1:L2) = PTEM2                                              
            END IF                                                              
            GO TO 810                                                           
         END IF                                                                 
  690    CONTINUE                                                               
C                                                                               
         DO 740 J=1,NSPEC                                                       
         XJ =INDXS(J)                                                           
         IF (RE(M1,NR).EQ.SPECI(J)) THEN                                        
            PTEM3 = '+K('//XR//')*Y('//XN//')*Y('//XJ//')*D'                    
            L2 = LTMP3+LTP(NS)                                                  
            LTP(NS) =  L2                                                       
            L1 = LP(NS)+1                                                       
            PROD(L1:L2) = PTEM3                                                 
            GO TO 810                                                           
         END IF                                                                 
 740     CONTINUE                                                               
 750     CONTINUE                                                               
C                                                                               
         DO 800 JC=1,INUM1                                                      
         XC = INDXC(JC)                                                         
         DO 765 M=1,2                                                           
         IF (RE(M,NR).EQ.SPSET(JC)) THEN                                        
            IF (M.EQ.1) THEN                                                    
               M1=M+1                                                           
            ELSE                                                                
               M1=M-1                                                           
            END IF                                                              
            GO TO 785                                                           
         END IF                                                                 
 765     CONTINUE                                                               
         GO TO 800                                                              
 785     DO 790 KC=1,INUM1                                                      
         XC1 = INDXC(KC)                                                        
         IF (RE(M1,NR).EQ.SPSET(KC)) THEN                                       
            IF (RE(M1,NR).EQ.'PHOTON'.OR.RE(M1,NR).EQ.'CRP'
     *         .OR.SPSET(IC).EQ.'CRPHOT') THEN               
               PTEM4 = '+K('//XR//')*X('//XC//')'                               
               L2 = LTMP4+LTP(NS)                                               
               LTP(NS) = L2                                                     
               L1 = LP(NS)+1                                                    
               PROD(L1:L2) = PTEM4                                              
C           ELSE IF (RE(M1,NR).EQ.'   ') THEN                                   
C              PTEM6 = '+K('//XR//')*X('//XC//')*D'                             
C              L2 = LTMP6+LTP(NS)                                               
C              LTP(NS) = L2                                                     
C              L1 = LP(NS)+1                                                    
C              PROD(L1:L2) = PTEM6                                              
            ELSE                                                                
               PTEM5 = '+K('//XR//')*X('//XC//')*X('//XC1//')*D'                
               L2 = LTMP5+LTP(NS)                                               
               LTP(NS) = L2                                                     
               L1 = LP(NS)+1                                                    
               PROD(L1:L2) = PTEM5                                              
            END IF                                                              
            GO TO 810                                                           
         END IF                                                                 
 790     CONTINUE                                                               
 800     CONTINUE                                                               
 810     LP(NS) = LTP(NS)                                                       
      END IF                                                                    
 815  CONTINUE                                                                  
 820  CONTINUE                                                                  
C-----------------------------------------------------------------------        
C WRITE PRODUCTION TERMS TO OUTPUT FILE AS A CHECK ON FAILURE - THIS            
C WILL ONLY BE PERFORMED IF THE 'C' IN COLUMN 1 IS REPLACED WITH A BLANK        
C IN THE LINES UP TO STATEMENT NUMBER 870.                                      
C                                                                               
C     IF (LTP.GT.40) THEN                                                       
C     WRITE(12,840)NS,PROD(:40)                                                 
C 840 FORMAT(6X,'PROD(',I3,')= ',A40)                                           
C     NLINE = (LTP(NS)-40)/NPERL                                                
C     IF ((LTP(NS)-40).GT.NLINE*NPERL) NLINE = NLINE+1                          
C     DO 860 J=1,NLINE                                                          
C     J1 = J*40+(4*(J-1)+1)+(J-1)*10                                            
C     J2 = J1+NPERL-1                                                           
C     WRITE(12,1006)PROD(J1:J2)                                                 
C 860 CONTINUE                                                                  
C     ELSE                                                                      
C     WRITE(12,840)NS,PROD(:40)                                                 
C 870 END IF                                                                    
C-----------------------------------------------------------------------        
      IF (LTL(NS).EQ.0) LTL(NS) = 1                                             
      IF (LTP(NS).EQ.0) LTP(NS) = 1                                             
      L3 = LTL(NS)                                                              
      L4 = LTP(NS)                                                              
      IF (LTL(NS).EQ.1) THEN                                                    
         PROD(1:1) = ' '                                                        
         L5 = L4                                                                
      ELSE IF (LTP(NS).EQ.1) THEN                                               
         L5 = L3+13                                                             
         PROD(1:L5) = '+(Y('//XS//')*('//LOSS(1:L3)//'))'                       
      ELSE                                                                      
         L5 = L3+L4+13                                                          
         PROD(L4+1:L5) = '+(Y('//XS//')*('//LOSS(1:L3)//'))'                    
      END IF                                                                    
C                                                                               
C OUTPUT SECTION                                                                
C                                                                               
      L6 = L5+1                                                                 
      PROD(L6:L6+NPERL) = ' '                                                   
      IF (L5.LE.LCHAR1) THEN                                                    
         IF (L5.GT.40) THEN                                                     
            WRITE(12,1007)NS,PROD(:40)                                          
            NLINE = (L5-40)/NPERL                                               
            IF ((L5-40).GT.NLINE*NPERL) NLINE=NLINE+1                           
            DO 925 J = 1,NLINE                                                  
            J1 = J*40+(4*(J-1)+1)+(J-1)*10                                      
            J2 = J1+NPERL-1                                                     
            WRITE(12,1006)PROD(J1:J2)                                           
 925        CONTINUE                                                            
            GO TO 1000                                                          
         ELSE                                                                   
            WRITE(12,1007)NS,PROD(:40)                                          
            GO TO 1000                                                          
         END IF                                                                 
      ELSE                                                                      
C                                                                               
C THIS ELSE-END IF CLAUSE IS ONLY FOR SUBDIVISION OF EQUATIONS WHOSE            
C CONTINUATION LINES EXCEED 18 BEFORE OUTPUTTING                                
C                                                                               
         I1 = 1                                                                 
         I2 = LCHAR                                                             
         IGROUP = L5/LCHAR                                                      
         IF ((IGROUP*LCHAR).NE.L5)  IGROUP = IGROUP+1                           
         DO 990 I=1,IGROUP                                                      
         IF (I.GT.1) I2 = I1+LCHAR                                              
         LASTL = I2-NPERL+25                                                    
         DO 935 J=1,26                                                          
         J1 = LASTL+J                                                           
         IF (PROD(J1:J1+2).EQ.'*(-') THEN                                       
            J2 = J1-I1-IBCK                                                     
            WS(1:J2) = PROD(I1:J1-IBCK)                                         
            J1 = J1-IBCK                                                        
            GO TO 940                                                           
         END IF                                                                 
         IF (PROD(J1:J1).EQ.'+') THEN                                           
            J2 = J1-I1                                                          
            WS(1:J2) = PROD(I1:J1-1)                                            
            GO TO 940                                                           
         END IF                                                                 
         IF (PROD(I1:I1).EQ.'+'.AND.PROD(J1:J1).EQ.'-') THEN                    
            J2 = J1-I1+2                                                        
            WS(1:J2) = PROD(I1:J1-1)//'))'                                      
            GO TO 940                                                           
         END IF                                                                 
         IF (PROD(I1:I1).EQ.'-'.AND.PROD(J1:J1).EQ.'-') THEN                    
            J2 = J1-I1+10                                                       
            WS(1:J2) = '+Y('//XS//')*('//PROD(I1:J1-1)//')'                     
            GO TO 940                                                           
         END IF                                                                 
 935     CONTINUE                                                               
 940     J3 = J2+1                                                              
         WS(J3:) = ' '                                                          
         I1 = J1                                                                
         LASTG = I-1                                                            
         IF (I.GT.1) THEN                                                       
            WRITE(12,945) I,LASTG                                               
 945        FORMAT(6X,'YD(',I2,') = YD(',I2,')')                                
         ELSE                                                                   
            WRITE(12,950) I                                                     
 950        FORMAT(6X,'YD(',I2,') = ')                                          
         END IF                                                                 
         NLINE = J2/NPERL                                                       
         IF ((NLINE*NPERL).LT.J2) NLINE =NLINE+1                                
         DO 955 K =1,NLINE                                                      
         K1 =(K-1)*IPL+K                                                        
         K2 =K1+IPL                                                             
         WRITE(12,1006) WS(K1:K2)                                               
 955     CONTINUE                                                               
         IF ((L5-I2).LE.LCHAR) THEN                                             
            WRITE(12,970) NS,I                                                  
 970        FORMAT(6X,'YDOT(',I3,') = +YD(',I2,')')                             
            L8 =L5-J1+1                                                         
            NLINE = L8/NPERL                                                    
            IF ((NLINE*NPERL).LT.L8) NLINE =NLINE+1                             
            K3 = J1-1                                                           
            IF (PROD(J1:J1).EQ.'-') THEN                                        
               WRITE(12,980)NS                                                  
 980           FORMAT(5X,'*',8X,'+(Y(',I3,')*(')                                
            END IF                                                              
            DO 985 K =1,NLINE                                                   
            K1 = (K-1)*IPL+K                                                    
            K2 =K1+IPL                                                          
            WRITE(12,1006) PROD(K3+K1:K3+K2)                                    
 985        CONTINUE                                                            
            GO TO 1000                                                          
         END IF                                                                 
990      CONTINUE                                                               
      END IF                                                                    
1000  CONTINUE                                                                  
C     CLOSE(UNIT=10)                                                            
C     CLOSE(UNIT=11)                                                            
C     CLOSE(UNIT=12)                                                            
      STOP                                                                      
 1001 FORMAT(//)                                                                
 1002 FORMAT(5(3X,A3,2X,A7,:))                                                
 1003 FORMAT(5(1X,'Y(',A3,')=',A,1X,:))                                         
 1004 FORMAT(1X,A4,1X,4(1A7,1X),A3,1X,A3,1PE8.2,1X,0PF5.2,1X,
     *       F8.1,A9)                
 1005 FORMAT(6X,'X(',I2,') = TOTAL(',I2,')',A40)                                
 1006 FORMAT(5X,'*',8X,A54)                                                     
 1007 FORMAT(6X,'YDOT(',I3,') = ',5X,A40)                                       
      END                                                                       

