C MAIN PROGRAM OF SXLSQI
C
C  LAST REVISED DEC. 23, 2003
C
C  BY C. F. BAES, JR.
C
C
C MAIN PROGRAM
C
C  THE MAIN PROGRAM ESTABLISHES THE SIZE OF TWO LINEAR ARRAYS
C
C     S AND N
C
C   THAT WILL STORE ARRAYS OF ADJUSTABLE LENGTH DETERMINED BY
C   THE SIZE OF THE ASSUMED MODEL AND THE DATA SET
C
C 
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION S(50000),N(800)
      COMMON/LGSN/LGS,LGN
C  ASSIGN THE LENGTH OF ARRAY S DECLARED ABOVE
      LGS=50000
C  ASSIGN THE LENGTH OF ARRAY N DECLARED ABOVE 
      LGN=800
C  CALL SUBROUTINE THAT WILL LOCATE IN S AND N VARIOUS ARRAYS TO
C   BE USED IN THE CALCULATIONS
      CALL ARRAYS(S,N)
C
      PAUSE
C
      END
C
C
C
C  SUBROUTINE TO ESTABLISH THE LENGTH AND POSITION OF VARIOUS ARRAYS
C   IN THE LINEAR ARRAYS S AND N
C
C   PARAMETERS THAT DETERMINE ARRAY SIZES
C
C         NAN = NUMBER OF AQUEOUS ANIONS
C         NCA = NUMBER OF AQUEOUS CATIONS
C         NNU = NUMBER OF NEUTRAL AQUEOUS SPECIES
C         NO  = NUMBER OF DATA POINTS
C         NPA = NUMBER OF ABSORBENCY PARAMETERS PER SPECIES
C         NPH = DENOTES PRESENCE (1) OR ABSENCE (0) OF ONE HEAT
C                PARAMETER FOR SPECIES
C         NSA = NUMBER OF PRODUCT SPECIES IN THE AQUEOUS PHASE
C         NSO = NUMBER OF PRODUCT SPECIES IN THE NONAQUEOUS PHASE
C         NUM = FLAG INDICATING THAT EFFECTS OF UNSYMMETRICAL MIXING
C                OF IONS WILL (1) OR WILL NOT (0) BE INCLUDED IN THE
C                PITZER TREATMENT
C         NUP = NUMBER OF USER-DEFINED ADJUSTABLE PARAMETERS
C         NVP = NUMBER OF PARAMETERS BEING REFINED
C
C     ARRAY         CONTAINS                LENGTH      STARTING POSITION 
C
C     AC     NONAQ.SOLUTE ACTIVITY          2+NSO           S(JACS)
C             (X SCALE)
C     ALP    PITZER PARAM., CA.-AN.         2*NCA*NAN       S(JALP)
C     AM     ARRAY USED BY ORGLS            NVP*(NVP+3)/2   S(NLAM)  
C     AMEV   ARRAY USED BY ORGLS            NVP*(NVP+3)/2   S(NLAMV)
C     B0     PITZER PARAM., CA.-AN.         NCA*NAN         S(JB0)      
C     B1     PITZER PARAM., CA.-AN.         NCA*NAN         S(JB1)      
C     B2     PITZER PARAM., CA.-AN.         NCA*NAN         S(JB2)      
C     BP     PITZER PARAM., CA.-AN.         NCA*NAN         S(JBP)      
C     BG     PITZER PARAM., CA.-AN.         NCA*NAN         S(JBG)  
C     BGP    PITZER PARAM., CA.-AN.         NCA*NAN         S(JBGP)  
C     CAM    ANION CONCENTRATION            NAN             S(JAM)
C     CAP    CALC.AQ.PROD.SPC.CONC.         NSA             S(JCA)
C     CBM    NEUTRAL SP. CONCENTRATION      NNU             S(JBM)
C     CEOC   EXPON.IN NOM-AQ.COMB.RCTN.     8+NSO           S(JCEO)
C     CCM    CATION CONCENTRATION           NCA             S(JCM)
C     CG     PITZER PARAM., CA.-AN.         NCA*NAN         S(JCG)
C     COP    CALC.NONAQ.PROD.SPC.CONC.      NSO             S(JCO)
C     CP     PITZER PARAM., CA.-AN.         NCA*NAN         S(JCP)      
C     CS     EST.EQUIL.CONCENTRATIONS       8+NSO+NSA       S(JLCS)        
C     DCA    PITZER PARAM., NU.             NNU             S(JDCA)      
C     DIAG   ARRAY USED BY ORGLS            NVP             S(NLDG)
C     DP     ADJ.PARAMETER INCREMENTS       LP              S(JDP)
C     DV     ARRAY USED BY ORGLS            NVP             S(NLDV)
C     EAC    RCT.CON.EXP., AQ.EQUIL.CONST.  8*NSA           S(JEAC)
C     EOC    RCT.CON.EXP., NONAQ.EQ.CONST.  8*NSO           S(JEOC)
C     EP     ARRAY USED BY SMPLX            8+NSO+NSA       S(JLEP)
C     ETHA   PITZ.UNSYM.MXG.PARAMS, ANIONS  NUM*NAN*NAN     S(JEA)  
C     ETHAP  PITZ.UNSYM.MXG.PARAMS, ANIONS  NUM*NAN*NAN     S(JEAP)
C     ETHC   PITZ.UNSYM.MXG.PARAMS, CATIONS NUM*NCA*NCA     S(JEC)  
C     ETHCP  PITZ.UNSYM.MXG.PARAMS, CATIONS NUM*NCA*NCA     S(JECP)
C     EVAL   ARRAY USED BY ORGLS            NVP             S(NLEVA)
C     EVEC   ARRAY USED BY ORGLS            NVP*NVP         S(NLEVC)
C     EXC    ARRAY USED BY SMPLX            8+NSO+NSA       S(JLEX)
C     FAC    FACTOR IN AQ.EQUIL.CONST.      8*NSA           S(JFAC)
C     FAK    AQ.PROD.SP.FORMATION CONST.    NSA             S(JFAK)
C     FCL    FCTNS.OF XCL TO BE SATISFIED   8+NSO+NSA       S(JLFCL)
C     FOC    FACTOR IN NONAQ.EQ.CONST.      8*NSO           S(JFOC)
C     FOK    NONAQ.PROD.SP.FORM'N.CONST.    NSO             S(JFOK)
C     FX     ARRAY USED BY SMPLX            11+NSO+NSA      S(JLFX)
C     G      ACTIVITY COEFFICIENTS          8+NSO+NSA       S(JLG)
C     GO     SOLUBILITY PARAMETERS          2+NSO           S(JLGO)
C     GX     NONAQ.MOLE FRACT.ACT.COEF.     2+NSO           S(JGXS)
C     GX0    LIMITING ACTIVITY COEF.        2+NSO           S(JGX0)
C     IAN    ANION IDENTIFIER               NAN             N(IIAN)
C     ICA    CATION IDENTIFIER              NCA             N(IICA)
C     INU    NEUTRAL SP. INDETIFIER         NNU             N(IINU)
C     IOUT   ARRAY USED BY ORGLS            NVP             N(IOU)
C     KI     REFINEMENT INTEGERS            LP              N(IKI)
C     KX     ARRAY USED BY SMPLX            11+NSO+NSA      N(IKX)
C     NCAP   FMLA.COEFS., AQU.PROD.SPECIES  7*NSA           N(ICA)
C     NCOP   FMLA.COEFS., NONAQ.PROD.SP.    8*NSO           N(ICO)
C     NF     INDEX OF CONCNS.TO BE DETMD.   8+NSO+NSA       N(INF)
C     NPSI   NO.OF EST.CONCNS.              NSO             N(INPS)
C     NIDT   INTEGER DEFINING DATA TYPE     NO              N(IDT)
C     NZA    CHARGES ON AQ.SPECIES          6+NSA           N(IZA)
C     NZPO   CHARGE ON NONAQ.PROD.SPECIES   NSO             N(IZO)
C     P      ADJUSTABLE PARAMETERS          LP              S(JLP)
C     PAA    PITZER PARAM., AN.-AN.         NAN*NAN         S(JPAA)      
C     PB     ARRAY USED BY SMPLX            8+NSO+NSA       S(JLPB)
C     PCAA   PITZER PARAM., CA.-AN.-AN.     NCA*NAN*NAN     S(JPCAA)      
C     PCC    PITZER PARAM., CA.-CA.         NCA*NCA         S(JPCC)      
C     PCCA   PITZER PARAM., CA.-CA.-AN.     NCA*NCA*NAN     S(JPCCA)      
C     PD     ARRAY USED BY ORGLS            NVP             S(NLPD)
C     PH     VOLUME FRACTION                2+NSO           S(JPH)
C     PHV    ARRAY USED BY ORGAC            2+NSO           S(JPHV)
C     PNA    PITZER PARAM., NU.-AN.         NNU*NAN         S(JPNA)      
C     PNC    PITZER PARAM., NU.-CA.         NNU*NCA         S(JPNC)      
C     PNN    PITZER PARAM., NU.-NU.         NNU*NNU         S(JPNN)
C     PT     ARRAY USED BY SMPLX            (8+NSO+NSA)*    S(JLPT)
C                                           (11+NSO+NSA)  
C     RAC    FMLA.COEFS., AQU.PROD.SPECIES  8*NSA           S(JRAC)    
C     REOC   EXPON.IN NOM-AQ.COMB.RCTN.     (8+NSO)*NSO*    S(JREO)
C                                           (MSO-1)*5/2                                                            
C     RFOK   EQ.CONST.,NOM-AQ.COMB.RCTN.    NSO*(MSO-1)*5/2 S(JRFO)
C     ROC    FMLA.COEFS., NONAQ.PROD.SP.    8*NSO           S(JROC)
C     ROW    ARRAY USED BY ORGLS            NVP             S(NLRW)
C     RZIR   RATIO OF CHARGE ON NONAQ.ION   NSO             S(JRZI)
C             TO THAT ON REF.ION
C     SAM    ARRAY USED BY ORGLS            NVP             S(NLSM)
C     SCPSI   SUM OF ION CONCNS.            NSO             S(JSCP)
C     SEOC   EXPON.IN NOM-AQ.COMB.RCTN.     (8+NSO)*NSO*    S(JSEO)
C                                           (MSO-1)/2                                                            
C     SFOK   EQ.CONST.,NOM-AQ.COMB.RCTN.    NSO*(MSO-1)/2   S(JSFO)
C     SIGYO  UNCERTAINTY IN OBS.QUANT.      NO              S(JLS)
C     SV     MASSON COEF.SV OF AQU.IONS     6+NSA           S(JSV)
C     V      MOLAR VOL. OF SPECIES          8+NSO           S(JLMV)
C     VC     ARRAY USED BY ORGLS            NVP             S(NLV)
C     V0     MASSON COEF.V0 OF AQU.IONS     6+NSA           S(JVO)  
C     VRT    ARRAY USED BY ORGAC            2+NSO           S(JVRT)
C     W      M.W. OF SPECIES                8+NSO+NSA       S(JLMW)
C     WA     ARRAY USED BY HYBRD1           (8+NSO+NSA)*    S(JLWA)
C                                           (3*(8+NSO+NSA)+7)/2
C     X      INDEPENDENT VARIABLES          10*NO           S(JLX)
C     XC     LG.EQUIL.CONCS., CURR.DATA PT  8+NSO+NSA       S(JLXC)
C     XCL    LG.CURR.CONCS.TO BE DETMD.     8+NSO+NSA       S(JLXCL)
C     XCAS   LG.EQ.CONS., DATA PTS.W/O      NO*(6+NSA)      S(JLXCAS)
C             NONAQ.PHASE
C     XCOS   LG.EQ.CONS., DATA PTS.W/O      NO*(8+NSO+NSA)  S(JLXCOS)
C             AQ.PHASE
C     XCS    LG.EQ.CONS., W BOTH PHASES     NO*(8+NSO+NSA)  S(JLXCS)
C     XF     NONAQ.SOLUTE MOLE FRACTION     2+NSO           S(JXF)
C     XS     LG.GVN.CONCENTRATIONS.         11*NO           S(JLXS)
C     YO     OBSERVED QUANTITY              NO              S(JLY)
C     ZA     ANION CHARGE                   NAN             S(JZA)
C     ZC     CATION CHARGE                  NCA             S(JZC)
C     ZFA    ANION CHARGE FUNCTION          NAN             S(JZFA)
C     ZFC    CATION CHARGE FUNCTION         NCA             S(JZFC)
C
C
C  OTHER QUANTITIES:
C
C     SMINS = A POSITIVE VALUE THAT ESTABLISHES THE PROCEDURE TO BE
C              USED IN SOLVING SETS OF NON-LINEAR SIMULTANEOUS EQUATIONS.
C              WHEN THE RESIDUALS OF THESE EQUATIONS EXCEED A LIMIT
C              DERIVED FROM THE VALUE OF SMINS, A SIMPLEX PROCEDURE IS
C              USED TO REDUCE THE RESIDUALS BELOW THIS LIMIT. THEREAFTER,
C              A MORE EFFICIENT, BUT LESS ROBUST, EQUATION SOLVER (HYBRD2)
C              IS EMPLOYED. A VALUE OF SMINS = 2 IS RECOMMENDED. IF PROBLEMS
C              ARE ENCOUNTERED WITH CONVERGENCE IN SUBROUTINE HYDRD1 THAT
C              CAN'T BE OVERCOME ANY OTHER WAY, THE VALUE OF SMINS CAN BE
C              REDUCED, THOUGH THIS WILL INCREASE THE COMPUTATION TIME.
C
C     MXHCY = IN HYBRD2, THE MAXIMUM NUMBER OF ITERATIONS PERMITTED BEFORE
C              A CALCULATION IS TERMINATED. THE RECOMMENDED VALUE IS 1000
C
C      NTRY = IN HYBRD2, THE NUMBER OF ITERATIONS AFTER WHICH A CALCULATION
C              IS TERMINATED IF THE RESIDUAL OF THE SET OF EQUATIONS TO BE
C              SATISFIED IS NOT REDUCED BY BY A LEAST ONE PER CENT. THE VALUE
C              OF THIS QUANTITY DEPENDS ON THE NUMBER AND PROPERTIES OF THE
C              EQUATIONS BEING SOLVED. IT CAN BE AS SMALL AS 20 OR AS LARGE
C              FINDS NECESSARY. 
C
C        LP = NSO*(1+NPH+NPA)+2*(NPH+NPA)+NSA+2+NSO+4*NCA*NAN+
C              (1+NAN)*NCA*(NCA-1)/2+(1+NCA)*NAN*(NAN-1)/2+NNU+1+NNU*NCA+
C              NNU*NAN+NNU*(NNU+1)/2+NUP
C
C
C
C  ORDER IN WHICH QUANTITIES ARE STORED IN VARIOUS ARRAYS
C
C        ORDER                 ARRAY                   QUANTITY
C
C   M,N,A,B,X,Y,H/OH:   NCAP(1,J) - NCAP(7,J)  FORMULA COEF.,AQ.SP. J
C
C   M,N,A,B,X,Y,H/OH,   NCOP(1,J) - NCOP(8,J)  FORMULA COEF.,NONAQ.
c    H20:                                       SP. J
C
C   HA,B,M,N,X,Y,H/OH:  X(1,J) - X(7,J)        GIVEN CONC., POINT J
C
C   HA,B,M,N,X,Y,H/OH:  XS(1,J) - XS(7,J)      LOG GIVN.CONC., POINT J
C   M,N,X,Y:            XS(8,J) - XS(11,J)     SHIFTED VALS.IF NO AQ.PH.
C
C   HA,B,M,N,X,Y,H,OH,  XC(1) - XC(8),         LOG.EQ.SPEC.CONC. FOR
C    NONAQ.SP.1 - NSO,  XC(9) - XC(8+NSO)       A SINGLE DATA POINT
C    AQ.SP.1 - NSA:     XC(9+NSO) - XC(8+NSO+NSA)
C
C       (SAME)          XCS(1,J) -             LOG.EQ.SPEC.CONC. FOR
C                       XCS(8+NSO+NSA,J)        POINT J, WITH BOTH PHASES
C
C       (SAME)          XCOS(1,J) -            LOG.EQ.SPEC.CONC. FOR
C                       XCOS(8+NSO+NSA,J)       POINT J, W/O AQ.PHASE
C
C       (SAME)          CS(1) - CS(8+NSO+NSA)  EST.CONC.OF A SPECIES
C
C       (SAME)          G(1) - G(8+NSO+NSA)    ACT.COEF. FOR A SPECIES
C
C       (SAME)          W(1) - W(8+NSO+NSA)    MOL.WT. OF SPECIES
C
C       (SAME)          NF(1) - NF(NEQ)        INDEX NO.ASSGND.EACH.OF NEQ
C                                               SPC.CONCS.TO BE DTRMD.                                              
C
C   HA,B,M,N,X,Y,H,OH,  V(1) - V(8+NSO)        MOL.VOL. OF SPECIES
C    NAQ..SP.1 - NSO:
C
C       (SAME)          CEOC(1) - CEOC(8+NSO)  NO. OF REACTANT SPCS.IN
C                                               A FORM.RCTN.
C
C       (SAME)          REOC(1,J,K) -          NO. OF SPCS.IN COMB'N.
C                         REOC(8+NSO,J,K)       RCTN.J IN NONAQ.PHASE
C                                               THAT EXCLUDES AQ.REACT.
C                                               K.
C                                          
C       (SAME)          SEOC(1,J) -            NO. OF SPCS.IN COMB'N.
C                         REOC(8+NSO,J)         RCTN.J IN NONAQ.PHASE
C                                               THAT EXCLUDES ALL AQ.
C                                               REACTANT SPECIES
C                                          
C   HA,B,M,N,X,Y        TCS(1) - TCS(6)        GVN.CONC.OF REACTANT SPEC.
C
C   HA,B,M,N,X,Y,H,H2O  EAC(1,J) - EAC(8,J)    NO. OF REACTANT SPCS.IN
C                                               FORM.RCTN.OF AQ..PRD.SP.J
C
C       (SAME)          FAC(1,J) - FAC(8,J)    FACTOR OF REACTANT.SPCS.IN
C                                               FORM.CNST.OF AQ..PRD.SP.J
C
C       (SAME)          EOC(1,J) - EOC(8,J)    NO. OF REACTANT SPCS.IN
C                                               FORM.RCTN.OF NONAQ.PRD.
C                                               SP.J
C
C       (SAME)          FOC(1,J) - FOC(8,J)    FACTOR OF REACTANT.SPCS.
C                                               IN FORM.CNST.OF NONAQ.
C                                               PRD.SP.J
C                                              
C       (SAME)          RFOC(1,J) - RFOC(8,J)  FACTOR OF REACTANT.SPCS.
C                                               IN FORM.CNST.OF NONAQ.
C                                               PRD.SP.J
C                                              
C       (SAME)          ROC(1,J) - ROC(8,J)    FORMULA COEF., NONAQ.
C                                               SPEC.J
C
C   HA,B,M,N,X,Y,H      RAC(1,J) - RAC(7,J)    FORMULA COEF., AQ.SPEC.J
C                                             
C   M,N,X,Y,H,OH,       V0(1) - V0(6+NSA)      MASSON COEFFICIENTS
C    AQ.PD.SP.1 - NSA:  SV(1) - SV(6+NSA)       FOR AQUEOUS SPECIES
C
C       (SAME)          XCAS(1,J)-             LOG.EQ.SPEC.CONC. FOR
C                       XCAS(6+NSA,J)           POINT J, W/O NONAQ.PHASE     
C
C       (SAME)          CAP(1) - CAP(NSA)      CALC.EQ.SPEC.COMC.
C
C   HA,B                NCR(1), NCR(2)         DEG.OF ASSOCIATION
C       (SAME)          NCR(3), NCR(4)         NO.OF H2O MOL./AGGREGATE
C
C   M,N,H,X,Y,OH:       NZA(1) - NZA(6)        ION CHARGE. W/O SIGN
C   AQ.PD.SP.1 - NSA:   NZA(7) - NZA(6+NSA)    SP. CHARGE WITH SIGN
C
C   NONAQ.PD.SP.        NZPO(1) - NZPO(NSO)    SP. CHARGE WITH SIGN
C    1 - NSO:
C
C   HA,B,NONAQ.PD.SP.   GO(1) - GO(2+NSO)      SOLUBILITY PARAMETERS
C    1 - NSO:
C
C       (SAME)          COP(1) - COP(NSO)      CALC.EQ.SPEC.COMC.
C 
C       (SANE)          GX0(1) - GX0(2+NSO)    LIMITING ACTIVITY COEF.
C                       PH(1) - PH(2+NSO)       VOLUME FRACTION
C                       AC(1) - AC(2+NSO)       ACTIVITY (X SXALE)
C                       GX(1) - GX(2+NSO)       MOLE FRACTION ACT.COEF.
C                       XF(1) - XF(2+NSO)       MOLE FRACTION
C
C
      SUBROUTINE ARRAYS(S,N)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      CHARACTER *60 LTITLE
      EXTERNAL NSIZE
      COMMON/MCHPRS/EPS
      COMMON/RSDLS/SMINS
      COMMON/CYCLES/MXHCY,NTRY
      COMMON/LGSN/LGS,LGN
      COMMON/LGP/LP
      COMMON/MD/NSO,NCBO,IREF,NPH,NPA,NSA,NIO,NIA,NCA,NAN,NNU,NUP,NUM
      COMMON/TIT/LTITLE
      COMMON/DT/TC,NO,NOF,ICNR      
      COMMON/GRP/AFT,BFT,QFT,NCR(4),WS,VS,DS,GOS,WR,VR,GOR,RT
      COMMON/FIT/SIGMA,NC,NVP,LAM,LEV
      COMMON/WTR/P0,DW,DC,AP,RLKW      
      COMMON/EQ/DEN,AWG,AW,TWSA,TVSA,VWPA,VWPO,RMBA,RMBO
      COMMON/PRP/XFS,GXS,AS,RMCO,THSV,THS,TSAI,TSOI,TWSO,TVSO
      COMMON/ENTS/JENT
      DIMENSION S(LGS),N(LGN)
      DATA RJPC,R,TC0,LM1,LM7,LM8/
     &4.184D+00,8.31451D+0,273.15D+0,1,7,8/
      DATA AVN,BK,E,PI/
     &6.022045D+23,1.380662D-16,4.803242D-10,3.1415926D+0/
      DATA ZERO,ONE,TWO,THREE,FOUR,FIVE,EIGHT,TEN,THOU/
     &0.0D+0,1.0D+0,2.0D+0,3.0D+0,4.0D+0,5.0D+0,8.0D+0,1.0D+1,1.0D+3/
C
C ESTABLISH THE RELATIVE MACHINE PRECISION
C
      EPTST=ONE
1     EPSV=EPTST
      EPTST=EPTST/TWO
      IF(EPTST.GT.ZERO) GOTO 1
      EPS=DSQRT(EPSV)
C
      WRITE(9,110) EPS
110   FORMAT(/1X,'THE NACHINE PRECISION IS:',1PE12.3/)
C
      OPEN (UNIT = 10, FILE = 'SXLSQI.PAR',STATUS='OLD')
      OPEN (UNIT = 11, FILE = 'SXLSQI.DAT',STATUS='OLD')
C
C  READ TITLE, TEMPERATURE, AND NUMBER OF DATA POINTS
      READ(11,*) LTITLE,TC,NO
      OPEN (UNIT = 12, FILE = 'SXLSQI.LSQ',STATUS='UNKNOWN')
      OPEN (UNIT = 13, FILE = 'SXLSQI.FIT',STATUS='UNKNOWN')
      WRITE(13,1013) LTITLE
1013  FORMAT(1X,A60//
     &1X,' PT         YO         YC      YO-YC',
     &'      Sigma (YO-YC)/Sg       X(1)       X(2)'/
     &1X,'	              X(3)       X(4)',
     &'       X(5)       X(6)       X(7)       X(8)'/)
      OPEN (UNIT = 14, FILE = 'SXLSQI.SPC',STATUS='UNKNOWN')
      WRITE(14,1014) LTITLE
1014  FORMAT(1X,A60//
     &1X,'Pt.      C(HA/NA)       C(B/NB)          C(M)',
     &'          C(N)          C(x)'/
     &1X,'             C(Y)          C(H)         C(OH)'/
     &1X,'        C(naq.sp) ...  C(aq.sp)....'/)
      OPEN (UNIT = 15, FILE = 'SXLSQI.GMA',STATUS='UNKNOWN')
      WRITE(15,1015) LTITLE
1015  FORMAT(1X,A60//
     &1X,'Pt. g(HA/NA)  g(B/NB)     g(M)     g(N)',
     &'     g(X)     g(Y)     g(H)    g(OH)'/
     &1X,'     g(masp)...g(asp)....g(el)'/
     &1X,'    a(HA/NA)  a(B/NB)  a(nasp)... ',
     &'     (These activities refer to the pure'/   
     &1X,'      a(H2O)   a(slv)   a(slv,el) ',
     &'                  liquid standard state)'/)
      OPEN (UNIT = 16, FILE = 'SXLSQI.MBL',STATUS='UNKNOWN')
      WRITE(16,1016) LTITLE
1016  FORMAT(1X,A60//
     &1X,'IPT ICN     Den (aq)       TSAI      TVWSA',
     &'       VWPA         RMBA'/
     &1X,'                           TSOI      TVWSO       VWPO',
     &'         RMBO'/)
      OPEN (UNIT = 17, FILE = 'SXLSQI.ITR',STATUS='UNKNOWN')
      OPEN (UNIT = 18, FILE = 'SXLSQI.REF',STATUS='UNKNOWN')
      OPEN (UNIT = 19, FILE = 'SXLSQI.AUX',STATUS='UNKNOWN')
C  READ MODEL PARAMETERS THAT DETERMINE ARRAY SIZES, CONVERGENCE CRITERIA,
C   NO. OF CYCLES
      READ(10,*) NSO,NPH,NPA,NSA,NIO,NIA,NUM,NCA,NAN,NNU,NUP,
     &SMINS,MXHCY,NTRY,NC
C  EVALUATE THE PRODUCT R*T (CAL/MOL)
      T=TC0+TC
      RT=R*T/RJPC
C  CALCULATE NEEDED PROPERTIES OF WATER
      CALL WATER(TC,P0,DW,DC,RLKW)
C
C      WRITE(*,*) RLKW
C      PAUSE
C
C  CALCULATE THE DEBYE-HUCKEL CONSTANT
      AP=(ONE/THREE)*
     &DSQRT(TWO*PI*AVN*DW/THOU)*(E**2/(DC*BK*T))**(THREE/TWO)
      AFT=DSQRT(TWO*PI*AVN/THOU)*(E**2/BK)**(THREE/TWO)
      BFT=DSQRT(EIGHT*PI*AVN*E**2/(THOU*BK))
      QFT=E**2/(TWO*BK)
C
C      WRITE(9,210) AFT,BFT,QFT
C210  FORMAT(1X,'            AFT            BFT            QFT'/
C     &1X,1P3E15.8)
C
C  SET LENGTH OF THE P, DP, AND KI ARRAYS FOR ADJUSTABLE PARAMETERS
      LP=NSO*(1+NPH+NPA)+2*(NPH+NPA)+NSA+2+NSO+4*NCA*NAN+
     &(1+NAN)*(NCA*(NCA-1))/2+(1+NCA)*(NAN*(NAN-1))/2+NNU+NNU*NCA+
     &NNU*NAN+(NNU*(NNU+1))/2+NUP
C
C      WRITE(9,220) LP
C220   FORMAT(1X,'LENGTH OF P ARRAY IS',I4)
C
C  SET LENGTH OF ARRAYS FOR DEFINING SPECIES
      LZA=6+NSA
      LCO=8*NSO
      LCA=7*NSA
      LRO=8*NSO
      LRA=8*NSA
C  SET LENGTH OF ARRAYS DEFINING COMBINATION REACTIONS
      LCM=NSO*(NSO-1)/2
      LACM=5*LCM
      LSCM=(8+NSO)*LCM
      LASCM=5*LSCM
C  SET LENGTH OF ARRAYS INVOLVING AQUEOUS SPECIES
      LICA=NCA*NAN
      LICA2=2*LICA
      LICC=NCA*NCA
      LIEC=NUM*LICC
      LIAA=NAN*NAN
      LIEA=NUM*LIAA
      LINC=NNU*NCA
      LINA=NNU*NAN
      LINN=NNU*NNU
      LICCA=NCA*NCA*NAN
      LICAA=NCA*NAN*NAN
C  SET LENGTH OF STORAGE FOR DATA RELATED ARRAYS
      LX=10*NO
      LXS=11*NO
C  SET LENGTH OF VARIOUS OTHER ARRAYS
      LXV=8+NSO
      LXC=8+NSO+NSA
      LXCX=LXC+3
      LXX=LXC*LXCX
      LXCS=LXC*NO
      LXCOS=LXC*NO
      LXCAS=LZA*NO
      LGO=2+NSO
      LWA=LXC*(3*LXC+7)/2
C  ESTABLISH LOCATIONS IN INTEGER ARRAY N
      ISIZE=1
C   THE INTEGER ARRAYS DEFINING SPECIES
      IZA=NSIZE(ISIZE,LZA,LM7)
      IZO=NSIZE(ISIZE,NSO,LM1)
      ICO=NSIZE(ISIZE,LCO,LM8)
      ICA=NSIZE(ISIZE,LCA,LM7)      
C   THE INTEGER ARRAY DESIGNATING ADJUSTABLE PARAMETERS
      IKI=NSIZE(ISIZE,LP,LM1)
C   THE INTEGER ARRAY DEFINING DATA TYPE
      IDT=NSIZE(ISIZE,NO,LM1)
C   THE INTEGER ARRAYS STORING THE INDEXES OF SPECIE CONCNS.TO BE CALCD.
      INF=NSIZE(ISIZE,LXC,LM1)
      IICA=NSIZE(ISIZE,NCA,LM1)
      IIAN=NSIZE(ISIZE,NAN,LM1)
      IINU=NSIZE(ISIZE,NNU,LM1)
C   THE INTEGER ARRAY IN SMPLX
      IKX=NSIZE(ISIZE,LXCX,LM1)
C   ANY OTHER INTEGER ARRAYS
      INPS=NSIZE(ISIZE,NSO,LM1)
      LGNA=ISIZE-1
      WRITE(9,120) LGN,LGNA
120   FORMAT(//1X,'THE GIVEN LENGTH OF INTEGER ARRAY N IS',I6/
     &1X,'THE REQUIRED LENGTH IS NOW',I10//)
      IF(LGNA.GT.LGN) THEN
        PAUSE
        STOP
      END IF
C  ESTABLISH LOCATIONS IN LINEAR ARRAY S
      JSIZE=1
C   THE ADJUSTABLE PARAMETER ARRAY P
      JLP=NSIZE(JSIZE,LP,LM1)
C   PARAMETER INCREMENTS THAT WILL FORM ARRAY DP
      JDP=NSIZE(JSIZE,LP,LM1)
C   MOLECULAR WEIGHTS AND MOLAR VOLUMES OF SPECIES
      JLMW=NSIZE(JSIZE,LXC,LM1)
      JLMV=NSIZE(JSIZE,LXV,LM1)
C   THE MASSON PARAMETERS FOR AQUEOUS SPECIES
      JV0=NSIZE(JSIZE,LZA,LM1)
      JSV=NSIZE(JSIZE,LZA,LM1)
C   COEFFICIENTS RELATED TO PRODUCT SPECIES
      JROC=NSIZE(JSIZE,LRO,LM8)
      JRAC=NSIZE(JSIZE,LRA,LM8)
      JEOC=NSIZE(JSIZE,LRO,LM8)
      JEAC=NSIZE(JSIZE,LRA,LM8)
      JFOC=NSIZE(JSIZE,LRO,LM8)
      JFAC=NSIZE(JSIZE,LRA,LM8)
      JFOK=NSIZE(JSIZE,NSO,LM1)
      JFAK=NSIZE(JSIZE,NSA,LM1)
C   CONCENTRATIONS, ACTIVITY COEFFICIENTS, SOLUBILITY PARAMETERS
      JCA=NSIZE(JSIZE,NSA,LM1)
      JCO=NSIZE(JSIZE,NSO,LM1)
      JAM=NSIZE(JSIZE,NAN,LM1)
      JBM=NSIZE(JSIZE,NNU,LM1)
      JCM=NSIZE(JSIZE,NCA,LM1)
      JLG=NSIZE(JSIZE,LXC,LM1)
      JLGO=NSIZE(JSIZE,LGO,LM1)
      JGX0=NSIZE(JSIZE,LGO,LM1)
      JVRT=NSIZE(JSIZE,LGO,LM1)
      JPH=NSIZE(JSIZE,LGO,LM1)
      JPHV=NSIZE(JSIZE,LGO,LM1)
      JXF=NSIZE(JSIZE,LGO,LM1)
      JGXS=NSIZE(JSIZE,LGO,LM1)
      JACS=NSIZE(JSIZE,LGO,LM1)
C   CONSTANTS AND EXPONENTS RELATED TO COMBINATION REACTIONS
      JCEO=NSIZE(JSIZE,LRO,LM1)
      JSEO=NSIZE(JSIZE,LSCM,LM1)
      JSFO=NSIZE(JSIZE,LCM,LM1)
      JRFO=NSIZE(JSIZE,LACM,LM1)
      JREO=NSIZE(JSIZE,LASCM,LM1)
C   PITZER INTERACTION COEFS. AND REACTANT ARRAYS
      JB0=NSIZE(JSIZE,LICA,LM1)
      JB1=NSIZE(JSIZE,LICA,LM1)
      JB2=NSIZE(JSIZE,LICA,LM1)
      JBP=NSIZE(JSIZE,LICA,LM1)
      JBG=NSIZE(JSIZE,LICA,LM1)
      JBGP=NSIZE(JSIZE,LICA,LM1)
      JALP=NSIZE(JSIZE,LICA2,LM1)
      JCAM=NSIZE(JSIZE,NAN,LM1)
      JCBM=NSIZE(JSIZE,NNU,LM1)
      JCCM=NSIZE(JSIZE,NCA,LM1)
      JCAM=NSIZE(JSIZE,NAN,LM1)
      JCG=NSIZE(JSIZE,LICA,LM1)
      JCP=NSIZE(JSIZE,LICA,LM1)
      JPCC=NSIZE(JSIZE,LICC,LM1)
      JEC=NSIZE(JSIZE,LIEC,LM1)
      JECP=NSIZE(JSIZE,LIEC,LM1)
      JDCA=NSIZE(JSIZE,NNU,LM1)
      JPAA=NSIZE(JSIZE,LIAA,LM1)
      JEA=NSIZE(JSIZE,LIEA,LM1)
      JEAP=NSIZE(JSIZE,LIEA,LM1)
      JPCAA=NSIZE(JSIZE,LICAA,LM1)
      JPCCA=NSIZE(JSIZE,LICCA,LM1)
      JPNA=NSIZE(JSIZE,LINA,LM1)
      JPNC=NSIZE(JSIZE,LINC,LM1)
      JPNN=NSIZE(JSIZE,LINN,LM1)
      JZA=NSIZE(JSIZE,NAN,LM1)
      JZC=NSIZE(JSIZE,NCA,LM1)
      JZFA=NSIZE(JSIZE,NAN,LM1)
      JZFC=NSIZE(JSIZE,NCA,LM1)
C   DATA-RELATED VALUES
      JLY=NSIZE(JSIZE,NO,LM1)      
      JLS=NSIZE(JSIZE,NO,LM1)      
      JLX=NSIZE(JSIZE,LX,LM1)
      JLXS=NSIZE(JSIZE,LXS,LM1)      
C   STORAGE ARRAY FOR EQUILIBRIUM CONCENTRATIONS      
      JLXC=NSIZE(JSIZE,LXC,LM1)
      JLCS=NSIZE(JSIZE,LXC,LM1)      
      JLXCS=NSIZE(JSIZE,LXCS,LM1)      
      JLXCAS=NSIZE(JSIZE,LXCAS,LM1)      
      JLXCOS=NSIZE(JSIZE,LXCOS,LM1)      
      JLXCL=NSIZE(JSIZE,LXC,LM1)      
      JLFCL=NSIZE(JSIZE,LXC,LM1)
C   ARRAYS USED IN CALCULATION OF CONCENTRATIONS
      JLEX=NSIZE(JSIZE,LXC,LM8)
      JLWA=NSIZE(JSIZE,LWA,LM1)
      JLEP=NSIZE(JSIZE,LXC,LM1)
      JLFX=NSIZE(JSIZE,LXCX,LM1)
      JLPB=NSIZE(JSIZE,LXC,LM1)
      JLPT=NSIZE(JSIZE,LXX,LM1)
      JRZI=NSIZE(JSIZE,NSO,LM1)
      JSCP=NSIZE(JSIZE,NSO,LM1)
      LGSA=JSIZE-1
      WRITE(9,130) LGS,LGSA
130   FORMAT(//1X,'THE GIVEN LENGTH OF PARAMETER ARRAY S IS',I6/
     &1X,'THE REQUIRED LENGTH IS NOW',I6//)
      IF(LGSA.GT.LGS) THEN
        PAUSE
        STOP
      END IF
C  READ MODEL PARAMETERS
      CALL READPR
     &(S(JLP),S(JDP),S(JALP),S(JLMW),S(JLMV),S(JV0),S(JSV),      
     &N(IZA),N(IZO),N(ICO),N(ICA),N(IKI))
C  READ DATA
      CALL READAT
     &(S(JLMW),S(JLMV),S(JV0),S(JSV),S(JLY),S(JLS),S(JLX),
     &N(IDT),N(IZA),N(ICA))
C   CALCULATE M.W. AND M.V. OF SPECIES
      CALL CAMWVL
     &(S(JLMW),S(JLMV),N(IZA),N(ICO),N(ICA))
C   COUNT PARAMETERS TO BE VARIED
      NVP=0
      KMX=0
      J=IKI-1
      DO 10 I=1,LP
      J=J+1
      IF(N(J).EQ.0) GOTO 10
      IF(N(J).GT.KMX) KMX=N(J)
      NVP=NVP+1
10    CONTINUE
C   SET ALL NON-ZERO KI VALUES TO THE MAXIMUM KI VALUE
      J=IKI-1
      DO 20 I=1,LP
      J=J+1
      IF (N(J).EQ.0) GOTO 20
      N(J)=KMX
20    CONTINUE
C  SET LENGTH OF ARRAYS THAT ARE USED BY ORGLS
      LAM=NVP*(NVP+3)/2
      LEV=NVP*NVP
C  LOCATE ORGLS ARRAYS TO BE STORED IN ARRAY S
      NLAM=NSIZE(JSIZE,LAM,LM1)
      NLV=NSIZE(JSIZE,NVP,LM1)
      NLDV=NSIZE(JSIZE,NVP,LM1)
      NLDG=NSIZE(JSIZE,NVP,LM1)
      NLPD=NSIZE(JSIZE,NVP,LM1)
      NLRW=NSIZE(JSIZE,NVP,LM1)
      NLEVA=NSIZE(JSIZE,NVP,LM1)
      NLEVC=NSIZE(JSIZE,LEV,LM1)
      NLAMV=NSIZE(JSIZE,LAM,LM1)
      NLSM=NSIZE(JSIZE,NVP,LM1)
      JSIZE=JSIZE-1     
      IOU=NSIZE(ISIZE,NVP,LM1)
      ISIZE=ISIZE-1
      WRITE(9,120) LGN,ISIZE
      IF(ISIZE.GT.LGN) THEN
        PAUSE
        STOP
      END IF
      WRITE(9,130) LGS,JSIZE
      IF(JSIZE.GT.LGS) THEN
        PAUSE
        STOP
      END IF
C
C  MAKE ANY NEEDED INITIAL CALLS TO OTHER SUBROUTINES
C
      JENT=0
      CALL EST
     &(S(JLX),S(JLXS),S(JLXCS),S(JLXCOS),S(JLXCAS),S(JLXC),IPTX,
     &IDCX,ICONVX,NCALX,ISTX,IEQX,IEHX)      
      JENT=1
C
C  BEGIN CALCULATIONS
C
      CALL ORGLS
     &(S(JLP),S(JDP),S(NLV),S(NLDV),S(NLDG),S(NLPD),S(NLRW),
     &S(NLEVA),S(NLEVC),S(NLAMV),S(NLSM),S(JLMW),S(JLMV),S(JV0),
     &S(JSV),S(JLCS),S(JLG),S(JLGO),S(JGX0),S(JVRT),S(JPH),
     &S(JPHV),S(JXF),S(JGXS),S(JACS),S(JLY),S(JLS),S(JLX),
     &S(JLXS),S(JLXCS),S(JLXCOS),S(JLXCAS),S(JLXC),S(JROC),S(JRAC),
     &S(JEOC),S(JCEO),S(JREO),S(JSEO),S(JEAC),S(JFOC),S(JFAC),
     &S(JFOK),S(JRFO),S(JSFO),S(JFAK),S(JLXCL),S(JLFCL),S(JLEX),     
     &S(JLWA),S(NLAM),S(JB0),S(JB1),S(JB2),S(JCP),S(JALP),     
     &S(JBP),S(JBG),S(JBGP),S(JCG),S(JPCC),S(JEC),S(JECP),     
     &S(JPAA),S(JEA),S(JEAP),S(JPCCA),S(JPCAA),S(JDCA),S(JPNC),     
     &S(JPNA),S(JPNN),S(JCAM),S(JCBM),S(JCCM),S(JZA),S(JZFA),     
     &S(JZC),S(JZFC),S(JCA),S(JCO),S(JLEP),S(JLPT),S(JLFX),     
     &S(JLPB),S(JRZI),S(JSCP),N(IKX),N(IZA),N(IZO),N(ICO),     
     &N(ICA),N(IDT),N(INF),N(IICA),N(IIAN),N(IINU),N(IKI),
     &N(IOU),N(INPS))
      CLOSE (UNIT=10)
      CLOSE (UNIT=11)
      CLOSE (UNIT=12)
      CLOSE (UNIT=13)
      CLOSE (UNIT=14)
      CLOSE (UNIT=15)
      CLOSE (UNIT=16)
      CLOSE (UNIT=17)
      CLOSE (UNIT=18)
      CLOSE (UNIT=19)
      RETURN
      END
C
C
C
C  FUNCTION TO CALCULATE STARTING LOCATION OF EACH ARRAY
C
      FUNCTION NSIZE(ISIZE,LENGTH,LMIN)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      NSIZE=ISIZE
      IF(LENGTH.GT.0) ISIZE=ISIZE+LENGTH
      IF(LENGTH.LE.0) ISIZE=ISIZE+LMIN
      RETURN
      END
C
C
C
C  READ MODEL PARAMETERS AND REFINEMENT INTEGERS AND SET INCREMENTS
C
C
C                              GLOSSARY
C
C     IREF = INDEX OF CHARGED NONAQUEOUS PRODUCT SPECIES USED TO DEFINE
C             FORMATION REACTION FOR OTHER CHARGED SPECIES
C             IREF IS IN THE RANGE 0 TO NSO
C     NCBO = FLAG DENOTING NO. AND KIND OF NONAQ.CHARGED SPECIES:
C              0 = NONE
C             -1 = AT LEAST ONE ANION (NOT ACCEPTABLE)
C              1 = AT LEAST ONE CATION (NOT ACCEPTABLE)
C              2 = AT LEAST ONE CATION AND ONE ANION
C
C
      SUBROUTINE READPR(P,DP,ALP,W,V,V0,SV,NZA,NZPO,NCOP,NCAP,KI)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      COMMON/MD/NSO,NCBO,IREF,NPH,NPA,NSA,NIO,NIA,NCA,NAN,NNU,NUP,NUM
      COMMON/GRP/AFT,BFT,QFT,NCR(4),WS,VS,DS,GOS,WR,VR,GOR,RT
      COMMON/LGP/LP
      DIMENSION P(*),DP(*),ALP(2,*),W(*),V(*),V0(*),SV(*),KI(*),
     &NZA(*),NZPO(*),NCOP(8,*),NCAP(7,*)
      DATA DLK,DPH,DPA,DPV,DLG,DLP2,DLP3/
     &0.001D+0,0.01D+0,0.01D+0,0.01D+0,0.001D+0,0.001D+0,0.0001D+0/
      DATA ONE,ZERO,THOU/1.0D+0,0.0D+0,1.0D+3/
      WRITE(*,999)
999   FORMAT(/1X,'READING THE PARAMETER FILE')
      READ(10,*) (NCR(I),I=1,4)
      READ(10,*) (NZA(I),I=1,(6+NSA))
C
      K=0
      IF(NSO.EQ.0) GOTO 35
      NAP=1+NPH+NPA
C
C  READ ORGANIC SPECIE PARAMETERS
      NCBOF=0
      IREF=0
      DO 15 I=1,NSO
C   READ PARAMETERS AND REFINEMENT INTEGERS
      READ(10,*) (NCOP(J,I),J=1,8),NZPO(I),
     &(P(L),KI(L),L=(K+1),(K+NAP))
C
C   VERIFY CHARGE ASSIGNED TO NONAQUEOUS PRODUCT SPECIES
      NZPOC=NCOP(1,I)*NZA(1)+NCOP(2,I)*NZA(2)-NCOP(3,I)-
     &NCOP(5,I)*NZA(4)-NCOP(6,I)*NZA(5)+NCOP(7,I)
      IF(NZPOC.NE.NZPO(I)) THEN
        WRITE(9,1001) I
1001    FORMAT(/1X,'FOR NONAQUEOUS PRODUCT SPECIES ',I3,':'/
     &  3X,'CHARGE ASSIGNED IS NOT COMPATIBLE WITH FORMULA'/)
        PAUSE
	STOP
      END IF
C
      IF(NZPO(I).NE.0) THEN
C      THIS IS A CHARGED SPECIES
        IF(NCBO.EQ.0)THEN
C        THIS IS THE FIRST CHARGED SPECIES
          IF(NZPO(I).GT.0) NCBO=1
	  IF(NZPO(I).LT.0) NCBO=-1
        END IF
        IF(NCBO.EQ.1) THEN
C        A CATION WAS SPECIFIED PREVIOUSLY
          IF(NZPO(I).LT.0) NCBO=2
        END IF
        IF(NCBO.EQ.-1) THEN
C        AN ANION WAS SPECIFIED PREVIOUSLY
          IF(NZPO(I).GT.0) NCBO=2
        END IF
        IF(P(K+1).EQ.ZERO.AND.KI(K+1).EQ.0.AND.IREF.EQ.0) THEN
C        THIS IS THE REFERENCE ION
          IREF=I
C
C          WRITE(9,4321) IREF
C4321      FORMAT(/1X,'IREF =',I3)
C          PAUSE
C
        END IF
      END IF
C
C   SET PARAMETER INCREMENTS
      KT=K
      K=K+1
      DP(K)=DLK
      IF(NPH.EQ.0) GO TO 5
      K=K+1
      DP(K)=DPH
5     IF(NPA.EQ.0) GO TO 15
      DO 10 J=1,NPA
      K=K+1
      DP(K)=DPA
10    CONTINUE
15    CONTINUE
C
C   CHECK FOR REQUIRED CONDITIONS IF NONAQ. IONS PRESENT
      IF(NCBO.EQ.1) THEN
        WRITE(9,1002) 
1002    FORMAT(/1X,'IN THE NONAQUEOUS PHASE:'/
     &  3X,'SPECIES HAVE BEEN SPECIFIED WITH POSITIVE CHARGE'/
     &  3X,'BUT NONE HAVE BEEN SPECIFIED WITH NEGATIVE CHARGE'/)
        PAUSE
	STOP
      END IF
      IF(NCBO.EQ.-1) THEN
        WRITE(9,1003) 
1003    FORMAT(/1X,'IN THE NONAQUEOUS PHASE:'/
     &  3X,'SPECIES HAVE BEEN SPECIFIED WITH NEGATIVE CHARGE'/
     &  3X,'BUT NONE HAVE BEEN SPECIFIED WITH POSITIVE CHARGE'/)
        PAUSE
	STOP
      END IF
      IF(NCBO.EQ.2) THEN
C      CHECK FOR RFERENCE ION
        IF(IREF.EQ.0) THEN
          WRITE(9,1004) 
1004      FORMAT(/1X,'IN THE NONAQUEOUS PHASE,',
     &    ' CHARGED SPECIES HAVE BEEN SPECIFIED, '/
     &    3X,'BUT NONE HAVE BEEN SELECTED AS A REFERENCE ION',
     &    ' BY FIXING ITS LOK K AT ZERO'/)
        PAUSE
	STOP
	END IF
      END IF
C
      IF(NAP.EQ.1) GO TO 35
C   READ PARAMETERS FOR REACTANTS HA AND B
      DO 30 I=1,2
      READ(10,*) (P(L),KI(L),L=(K+1),(K+NAP-1))
      IF(NPH.EQ.0) GO TO 20
      KT=K
      K=K+1
      DP(K)=DPH
20    IF(NPA.EQ.0) GO TO 30
      DO 25 J=1,NPA
      K=K+1
      DP(K)=DPA
25    CONTINUE
30    CONTINUE
C
35    IF(NSA.EQ.0) GOTO 50
C  READ AQUEOUS SPECIE PARAMETERS
      DO 40 I=1,NSA
C   READ PARAMETERS AND REFINEMENT INTEGERS
      K=K+1
      READ(10,*) (NCAP(J,I),J=1,7),DMY,P(K),KI(K),V0(6+I),SV(6+I)
C
C    VERIFY CHARGE ASSIGNED TO AQUEOUS PRODUCT SPECIES
      NZPA=NCAP(1,I)*NZA(1)+NCAP(2,I)*NZA(2)-NCAP(3,I)
     &-NCAP(5,I)*NZA(4)-NCAP(6,I)*NZA(5)+NCAP(7,I)
      IF(NZPA.NE.NZA(6+I)) THEN
        WRITE(9,1005) I
1005    FORMAT(/1X,'FOR AQUEOUS PRODUCT SPECIES ',I3,':'/
     &  3X,'CHARGE ASSIGNED IS NOT COMPATIBLE WITH FORMULA'/)
        PAUSE
	STOP
      END IF 
C
C   SET PARAMETER INCREMENTS
      DP(K)=DLK
40    CONTINUE
C  READ NONAQUEOUS SOLUBILITY PARAMETERS
50    READ(10,*) GOS 
      DO 55 I=1,(NSO+2)
      K=K+1
      READ(10,*) P(K),KI(K) 
      DP(K)=DLG
55    CONTINUE
      M=0
C  READ AQUEOUS PITZER PARAMETERS
C   READ B0, B1, B2, CP, A1, AND A2
      DO 70 I=1,NCA
      DO 65 J=1,NAN
      M=M+1
      READ(10,*) (P(L),KI(L),L=(K+1),(K+4)),ALP(1,M),ALP(2,M)
      KT=K
      DO 60 L=1,4
      K=K+1
      DP(K)=DLP2
60    CONTINUE
65    CONTINUE
70    CONTINUE
C   READ THETA AND PSI
      IF(NCA.LT.2) GOTO 90
      DO 85 I=1,(NCA-1)
      DO 80 J=(I+1),NCA
      READ(10,*) (P(L),KI(L),L=(K+1),(K+2))
      KT=K
      K=K+1
      DP(K)=DLP3
      K=K+1
      DP(K)=DLP3
      IF(NAN.EQ.1) GOTO 80
      DO 75 L=1,(NAN-1)
      K=K+1
      READ(10,*) P(K),KI(K)
      DP(K)=DLP3
75    CONTINUE
80    CONTINUE
85    CONTINUE
90    IF(NAN.LT.2) GOTO 110
      DO 105 I=1,(NAN-1)
      DO 100 J=(I+1),NAN
      READ(10,*) (P(L),KI(L),L=(K+1),(K+2))
      KT=K
      K=K+1
      DP(K)=DLP3
      K=K+1
      DP(K)=DLP3
      IF(NCA.EQ.1) GOTO 100
      DO 95 L=1,(NCA-1)
      K=K+1
      READ(10,*) P(K),KI(K)
      DP(K)=DLP3
95    CONTINUE
100   CONTINUE
105   CONTINUE
110   IF(NNU.EQ.0) GOTO 140
C   READ D AND RHO
      DO 130 I=1,NNU
      READ(10,*) (P(L),KI(L),L=(K+1),(K+2))
      KT=K
      K=K+1
      DP(K)=DLP2
      K=K+1
      DP(K)=DLP2
      IF(NCA.LT.2) GOTO 120
      DO 115 J=2,NCA
      K=K+1
      READ(10,*) P(K),KI(K)
      DP(K)=DLP2
115   CONTINUE
120   DO 125 J=1,NAN
      K=K+1
      READ(10,*) P(K),KI(K)
      DP(K)=DLP2
125   CONTINUE
130   CONTINUE
C   READ LAMBDA VALUES
      DO 135 I=1,(NNU*(NNU+1)/2)
      K=K+1
      READ(10,*) P(K),KI(K)
      DP(K)=DLP2      
135   CONTINUE
140   IF(NUP.EQ.0) GOTO 150
C  READ USER-DEFINED PARAMETERS
      DO 145 I=1,NUP
      K=K+1
      READ(10,*) P(K),KI(K),DP(K)
145   CONTINUE
150   CONTINUE
C
      WRITE(9,1006) K,LP
1006  FORMAT(/1X,'AT THE END OF S. READPR,'/
     &3X,'THE NUMBER OF ADJUSTABLE PARAMETERS READ IS',I3/
     &3X,'THE MAXIMUM NUMBER PERMITTED IS', I3)
      IF(K.GT.LP) THEN
        PAUSE
	STOP
      END IF 
C
      RETURN
      END
C
C
C
C  READ THE DATA
C
C
      SUBROUTINE READAT(W,V,V0,SV,YO,SIGYO,X,NIDT,NZA,NCAP)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      DIMENSION W(*),V(*),V0(*),SV(*),YO(*),SIGYO(*),X(10,*),
     &NIDT(*),NZA(*),NCAP(7,*)
      COMMON/MD/NSO,NCBO,IREF,NPH,NPA,NSA,NIO,NIA,NCA,NAN,NNU,NUP,NUM
      COMMON/DT/TC,NO,NOF,ICNR      
      COMMON/GRP/AFT,BFT,QFT,NCR(4),WS,VS,DS,GOS,WR,VR,GOR,RT
      DATA IT1,IT2,IT3,IT4/10,100,1000,10000/
      DATA ONE,ZERO/1.0D+0,0.0D+0/
      DATA FMIN/1.0D-13/
      WRITE(*,910)
910   FORMAT(/1X,'READING THE DATA FILE')
      NOF=NO            
C  READ M.W., MOLAR VOL., AND DIELECTRIC CONSTANT OF NONAQ.SOLVENT;
C   M.W., MOLAR VOL., AND SOLUBILITY PARAMETER OF REFERENCE SOLUTE;
      READ(11,*) WS,VS,DS
      READ(11,*) WR,VR,GOR
C  READ M.W. AND NONAQUEOUS MOLAR VOLS.OF HA AND B
      READ(11,*) W(1),V(1)
      READ(11,*) W(2),V(2)
C  READ M.W., APP.MOL.VOLS., AND MASSON CONSTS. OF AQ.IONS
      DO 5 J=1,6
      READ(11,*) W(J+2),V(J+2),V0(J),SV(J)
5     CONTINUE
C
C  BEGIN READING THE DATA
C
      DO 25 J=1,NO
      READ(11,*) NIDT(J),YO(J),SIGYO(J),(X(I,J),I=1,10)
C   VERIFY THAT FIRST SIX VALUES OF X(I,J) ARE ZERO OR POSITIVE
      DO 10 I=1,6
      IF(X(I,J).LT.ZERO) THEN
        WRITE(9,1000) J,I 
1000  FORMAT(/
     &1X,'FOR DATA POINT',I3,':'/
     &3X,'REACTANT CONCENTRATION ',I1,' IS LESS THAN ZERO.'/)
        PAUSE
	STOP
      END IF
10    CONTINUE
C
C
C
C   CHECK THE PROPERTIES OF DATA POINT DEFINED BY FIRST ENTRY
C
C
C                IAO  FLAG INDICATING: (0) THE PRESENCE OF ONLY THE
C                     NONAQUEOUS PHASE; (1) THE PRESENCE OF BOTH PHASES; OR
C                     (2) THE PRESENCE OF ONLY THE AQUEOUS PHASE.
C
C                IEQ  FLAG INDICATING: (0) INITIAL TOTAL CONCNS. HAVE
C                     BEEN SUPPLIED FOR ALL REACTANTS OTHER THAN H+; OR
C		      (1) EQUILIBRIUM CONCNS. HAVE BEEN SUPPLIED FOR ALL
C                     REACTANTS OTHER THAN H+.
C
C                IEH  FLAG INDICATING: (0) INITIAL TOTAL ACID CONCN. HAS
C                     BEEN SUPPLIED; (1) THE EQUIL. ACID CONCN. HAS BEEN
C                     SUPPLIED; OR (>1) THE EQUILIBRIUM CONCENTRATION OF
C                     H+ HAS BEEN SUPPLIED.
C
C                ICN  FLAG INDICATING CONCENTRATION SCALE: (1) MOLARITY,
C                     (2) MOLALITY, OR (3) MOL/KG SOLN.
C
C
C    VERIFY THE MAGNITUDE THESE QUANTITIES
C
      NSVE=NIDT(J)
      IRES=NIDT(J)
      IY=IRES/IT4
      IRES=IRES-IT4*IY
      IAO=IRES/IT3
      IF(IAO.LT.0.OR.IAO.GT.2) THEN
        WRITE(9,1001) J 
1001  FORMAT(/
     &1X,'FOR DATA POINT',I3,':'/
     &3X,'INTEGER 3 OF FIRST ENTRY SHOULD BE 0, 1, OR, 2.'/)
        PAUSE
	STOP
      END IF
      IRES=IRES-IT3*IAO
      IEQ=IRES/IT2
      IF(IEQ.LT.0.OR.IEQ.GT.1) THEN
        WRITE(9,1002) J 
1002  FORMAT(/
     &1X,'FOR DATA POINT',I3,':'/
     &3X,'INTEGER 4 OF FIRST ENTRY SHOULD BE 0 OR 1.'/)
        PAUSE
	STOP
      END IF
      IRES=IRES-IT2*IEQ
      IEH=IRES/IT1
      IF(IEH.LT.0.OR.IEH.GT.5) THEN
        WRITE(9,1003) J 
1003  FORMAT(/
     &1X,'FOR DATA POINT',I3,':'/
     &3X,'INTEGER 5 OF FIRST ENTRY SHOULD BE 0, 1, 2, 3, 4, OR 5.'/)
        PAUSE
     	STOP
      END IF
      ICN=IRES-IT1*IEH
C    SAVE THE INITIAL VALUE OF ICN TO ESTABLISH THE CONCENTRATION SCALE
C     OF MODELS
      IF(J.EQ.1) ICNR=ICN
      IF(ICN.LT.1.OR.ICN.GT.3) THEN
        WRITE(9,1004) J 
1004  FORMAT(/
     &1X,'FOR DATA POINT',I3,':'/
     &3X,'INTEGER 6 OF FIRST ENTRY SHOULD BE 1, 2, OR 3.'/)
        PAUSE
	STOP
      END IF
C
C    VALIDATE VALUES ASSIGNED TO IAO, IEQ, AND IEH:
C
C     WHEN VALUES ARE ASSIGNED TO ONE OR TWO OF THESE, THE VALUE(S)
C      ASSIGNED TO THE REMAINING ONE(S) IS SOMETIMES SPECIFIED.
C      REQUIRED VALUES ARE SHOWN IN PARENTHESES IN THE FOLLOWING TABLE.
C
C
C	               IAO     IEQ     IEH
C	             ------------------------
C			0      (0)     (0)
C			1	1      (³1)
C			2      (0)     (0 OR >1)
C	             ------------------------
C
C     IF UNACCEPTABLE VALUES ARE FOUND, EITHER THEY WILL BE CORRECTED
C      OR EXECUTION OF THE PROGRAM WILL BE TERMINATED
C
      IF(IAO.EQ.0) THEN
        IF(X(3,J).GT.ZERO.OR.X(4,J).GT.ZERO.OR.X(5,J).GT.ZERO.OR.
     &  X(6,J).GT.ZERO.OR.X(7,J).NE.ZERO) THEN
          IF(X(8,J).LT.ONE) THEN
	    WRITE(9,1010) J 
1010  FORMAT(
     &1X,'FOR DATA POINT',I3,':'/
     &3X,'THE AQ.PHASE IS ABSENT AND AQUEOUS REACTANTS',
     &' ARE PRESENT;'/
     &3X,'HENCE, A WATER ACTIVITY OF UNITY SHOULD BE ASSIGNED')
            PAUSE
            STOP
          END IF
        END IF
      END IF
      IF(IY.EQ.40.AND.IAO.NE.0) THEN
        WRITE(9,1011) J 
1011  FORMAT(/
     &1X,'FOR DATA POINT',I3,':'/
     &3X,'THE AQ.PHASE IS ASSUMED ABSENT'/)
	PAUSE
        IAO=0
      END IF
      IF(IY.EQ.41.AND.IAO.NE.1) THEN
        WRITE(9,1012) J
1012  FORMAT(/
     &1X,'FOR DATA POINT',I3,':'/
     &3X,'BOTH PHASES ARE ASSUMED PRESENT'/)
	PAUSE
        IAO=1
      END IF
      IF(IAO.EQ.0) THEN
C      NO AQUEOUS PHASE IS PRESENT. AN INFINITESIMAL AMOUNT
C       OF PURE WATER WILL BE ASSUMED TO BE EQUILIBRATED WITH THE. 
C       NONAQUEOUS PHASE
        IEQ=0 
	IEH=0
	X(7,J)=ZERO
      END IF
      IF(IAO.EQ.1) THEN
C      BOTH PHASES ARE PRESENT
        IF(IEQ.EQ.1.AND.IEH.EQ.0) THEN
          WRITE(9,1014) J
1014  FORMAT(/
     &1X,'FOR POINT NO.',I3,':'/
     &3X,'WITH BOTH PHASES PRESENT, IF THE EQUILIBRIUM CONCNS. OF '/
     &1X,'THE OTHER AQ.COMPONENTS ARE SPECIFIED, THEN EITHER THE '/
     &1X,'EQUIL. CONCN. OF TOTAL ACID OR BASE, OR THE EQUIL. CONCN. '/
     &1X,'OF H+ MUST BE SPECIFIED.')
          PAUSE
          STOP
        END IF
      END IF
      IF(IAO.EQ.2) THEN
C      NO NONAQUEOUS PHASE IS PRESENT.
        IEQ=0 
        IF(IEH.EQ.1) THEN
          WRITE(9,1015) J
1015  FORMAT(/
     &1X,'FOR POINT NO.',I3,':'/
     &3X,'WITH ONLY THE AQ. PHASE PRESENT, THE INITIAL CONCNS. OF '/     									
     &1X,'THE OTHER AQ.COMPONENTS MUST BE SPECIFIED, AND EITHER THE '/
     &1X,'INITIAL CONCN. OF TOTAL ACID OR BASE, OR THE EQUIL. CONCN. '/
     &1X,'OF H+ MUST BE SPECIFIED.')
          PAUSE
          STOP
        END IF
      END IF
C  IF EQUILIBRIUM AQ. COMPN. IS SPECIFIED, AQ. PROD. SPECIES INVOLVING
C   A OR B ARE NOT PERMITTED
      IF(IEQ.EQ.1.AND.NSA.GT.0) THEN
        DO 15 I=1,NSA
        IF(NCAP(3,I).GT.0.OR.NCAP(4,I).GT.0) THEN
	  WRITE(9,1016) J
1016  FORMAT(/
     &1X,'FOR POINT NO.',I3,':'/
     &3X,'IF THE EQUILIBRIUM AQUEOUS COMPOSITION IS SPECIFIED, THEN '/
     &1X,'NO AQ. PRODUCT SPECIES INVOLVING A OR B IS PERMITTED')
          PAUSE
          STOP
        END IF
15      CONTINUE
      END IF
      NIDT(J)=ICN+IT1*IEH+IT2*IEQ+IAO*IT3+IT4*IY
      IF(NIDT(J).NE.NSVE) THEN
        IF(NSVE.LT.100000) WRITE(9,1017) J,NSVE,NIDT(J)
1017  FORMAT(/
     &1X,'FOR DATA POINT',I3,':'/
     &3X,'THE INITIAL ENTRY HAS BEEN CHANGED FROM 0',I5,' TO 0',I5/)
        IF(NSVE.GE.100000) WRITE(9,1018) J,NSVE,NIDT(J)
1018  FORMAT(/
     &1X,'FOR DATA POINT',I3,':'/
     &3X,'THE INITIAL ENTRY HAS BEEN CHANGED FROM ',I6,' TO ',I6/)
        PAUSE
      END IF
C
      IF(IEH.LE.1.AND.IAO.NE.0) THEN
C      ESTABLISH CHARGE BALANCE IN THE AQUEOUS PHASE BY ADJUSTING ANION
C       CONCENTRATIONS
        IF(NZA(4).EQ.0) X(5,J)=ZERO
        IF(NZA(5).EQ.0) X(6,J)=ZERO
        IF(NZA(6).EQ.0.AND.X(7,J).LT.ZERO) X(7,J)=ZERO
	XADJ=ONE
20      FCB=X(3,J)*NZA(1)+X(4,J)*NZA(2)-X(5,J)*NZA(4)-
     &  X(6,J)*NZA(5)+X(7,J)
        IF(DABS(FCB).GE.FMIN) THEN
          WRITE(9,1019) J,NZA(5),X(6,J),NZA(4),X(5,J),NZA(6),
     &    X(7,J),FCB
1019  FORMAT(/
     &1X,'FOR DATA POINT',I3,':'/
     &3X,'ION CHARGES DO NOT BALANCE.'/
     &1X,'  Z(Y)   C(Y)     Z(X)   C(X)    Z(OH)   C(H/OH)    ',
     &'CHG. RES.'/
     &1X,2(I4,F11.7,1X),I4,F12.7,3X,1P1E10.2/
     &1X,'[RTN] TO CONTINUE'/)
          PAUSE
          IF(NZA(5).NE.ZERO.AND.XADJ.NE.ZERO) THEN
C          Y IS PRESENT: ADJUST ITS TOTAL CONC.
            XADJ=X(6,J)+FCB/NZA(5)
            IF(XADJ.LE.ZERO) XADJ=ZERO 
            WRITE(9,1020) X(6,J),XADJ
1020  FORMAT(/1X,'ADJUSTING C(Y) FROM',F11.7,' TO',F11.7/
     &1X,'[RTN] TO CONTINUE'/)
            PAUSE
            X(6,J)=XADJ
            GOTO 20
          END IF
          IF(NZA(4).NE.ZERO.AND.XADJ.NE.ZERO) THEN
C          X IS PRESENT: ADJUST ITS TOTAL CONC.
            XADJ=X(5,J)+FCB/NZA(4)
            IF(XADJ.LE.ZERO) XADJ=ZERO 
      WRITE(9,1021) X(5,J),XADJ
1021  FORMAT(/1X,'ADJUSTING C(X) FROM',F11.7,' TO',F11.7/
     &1X,'[RTN] TO CONTINUE'/)
            PAUSE
            X(5,J)=XADJ
            GOTO 20
          END IF
          IF(X(7,J).LT.ZERO.AND.XADJ.NE.ZERO) THEN
C          ADJUST TOTAL CONC. OF OH-
            XADJ=X(7,J)-FCB
            IF(XADJ.GE.ZERO) XADJ=ZERO
            CNOH=-X(7,J)
            CNOHR=-XADJ
            WRITE(9,1022) CNOH,CNOHR
1022  FORMAT(/1X,'ADJUSTING C(OH) FROM',F11.7,' TO',F11.7/
     &1X,'[RTN] TO CONTINUE'/)
            PAUSE
            X(7,J)=XADJ
            GOTO 20
          END IF
          WRITE(9,1023) J
1023  FORMAT(/1X,'AQ.ION CHARGE CAN NOT BE BALANCED FOR POINT',
     &I4/1X,'[RTN] TO CONTINUE'/)
          PAUSE
          STOP
        END IF
      END IF
25    CONTINUE
C
C  END READING THE DATA
C
C  CALCULATE SOLUTION COMPOSITIONS FOR TITRATION DATA
      DO 35 J=1,NO
      NDTY=NIDT(J)/IT4
      IF(NDTY.NE.30) GOTO 35
      IF(X(9,J).EQ.ZERO) THEN
        JT=J
	NOF=NOF-1
      END IF
      IF(X(10,J).EQ.ZERO) THEN
        J0=J
	NOF=NOF-1
      END IF
      IF(J.GT.J0.AND.J.GT.JT) THEN
        DO 30 I=1,6
        X(I,J)=(X(I,J0)*X(9,J)+X(I,JT)*X(10,J))/(X(9,J)+X(10,J))    
30      CONTINUE
      END IF
35    CONTINUE
      RETURN
      END
C
C
C
C CALCULATE M.W. AND MOLAR VOLUMES OF PRODUCT SPECIES
C
      SUBROUTINE CAMWVL(W,V,NZA,NCOP,NCAP)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      COMMON/MD/NSO,NCBO,IREF,NPH,NPA,NSA,NIO,NIA,NCA,NAN,NNU,NUP,NUM
      COMMON/GRP/AFT,BFT,QFT,NCR(4),WS,VS,DS,GOS,WR,VR,GOR,RT
      DIMENSION W(*),V(*),NZA(*),NCOP(8,*),NCAP(7,*)
      DATA ONE,ZERO/1.0D+0,0.0D+0/
      IF(NSO.EQ.0) GOTO 20
C  NONAQUEOUS PRODUCT SPECIES
      DO 10 I=1,NSO
      K=8+I
      IF(NCOP(7,I).GE.0) VTM=(V(7)+V(8))*NCOP(8,I)+V(7)*NCOP(7,I)
      IF(NCOP(7,I).LT.0) VTM=(V(7)+V(8))*NCOP(8,I)-V(8)*NCOP(7,I) 
      V(K)=V(3)*NCOP(1,I)+V(4)*NCOP(2,I)+(V(1)-V(7))*NCOP(3,I)+
     &V(2)*NCOP(4,I)+V(5)*NCOP(5,I)+V(6)*NCOP(6,I)+VTM
      IF(NCOP(7,I).GE.0) WTM=(W(7)+W(8))*NCOP(8,I)+W(7)*NCOP(7,I)
      IF(NCOP(7,I).LT.0) WTM=(W(7)+W(8))*NCOP(8,I)-W(8)*NCOP(7,I) 
      W(K)=W(3)*NCOP(1,I)+W(4)*NCOP(2,I)+(W(1)-W(7))*NCOP(3,I)+
     &W(2)*NCOP(4,I)+W(5)*NCOP(5,I)+W(6)*NCOP(6,I)+WTM
10    CONTINUE
C  AQUEOUS PRODUCT SPECIES
20    IF(NSA.EQ.0) GOTO 40
      DO 30 I=1,NSA
      K=8+NSO+I
      IF(NCAP(7,I).GE.0) WTM=W(7)*NCAP(7,I)
      IF(NCAP(7,I).LT.0) WTM=-W(8)*NCAP(7,I)
      W(K)=W(3)*NCAP(1,I)+W(4)*NCAP(2,I)+(W(1)-W(7))*NCAP(3,I)+
     &W(2)*NCAP(4,I)+W(5)*NCAP(5,I)+W(6)*NCAP(6,I)+WTM
30    CONTINUE
C  NONAQUEOUS EXTRACTANT SPECIES
40    V(1)=V(1)*NCR(1)+(V(7)+V(8))*NCR(3)
      W(1)=W(1)*NCR(1)+(W(7)+W(8))*NCR(3)
      V(2)=V(2)*NCR(2)+(V(7)+V(8))*NCR(4)
      W(2)=W(2)*NCR(2)+(W(7)+W(8))*NCR(4)
      RETURN
      END
C
C
C
C  THIS SUBROUTINE IS A SLIGHTLY MODIFIED VERSION OF THE GENERAL
C  LEAST-SQUARES PROGRAM ORGLS BY W. R. BUSING AND AND H. A. LEVY
C
C
      SUBROUTINE ORGLS(P,DP,VC,DV,DIAG,PD,ROW,
     &EVAL,EVEC,AMEV,SAM,W,V,V0,SV,CS,G,GO,GX0,VRT,PH,
     &PHV,XF,GX,AC,YO,SIGYO,X,XS,XCS,XCOS,XCAS,XC,ROC,RAC,
     &EOC,CEOC,REOC,SEOC,EAC,FOC,FAC,FOK,RFOK,SFOK,FAK,XCL,FCL,EXC,
     &WA,AM,B0,B1,B2,CP,ALP,BP,BG,BGP,CG,PCC,ETHC,ETHCP,
     &PAA,ETHA,ETHAP,PCCA,PCAA,DCA,PNC,PNA,PNN,CAM,CBM,CCM,ZA,ZFA,
     &ZC,ZFC,CAP,COP,EP,PT,FX,PB,RZIR,SCPSI,KX,NZA,NZPO,NCOP,
     &NCAP,NIDT,NF,ICA,IAN,INU,KI,IOUT,NPSI)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE     
      CHARACTER*10 ISTOP
      COMMON/LGP/NP
      COMMON/MD/NSO,NCBO,IREF,NPH,NPA,NSA,NIO,NIA,NCA,NAN,NNU,NUP,NUM
      COMMON/TIT/LTITLE
      COMMON/DT/TC,NO,NOF,ICNR      
      COMMON/GRP/AFT,BFT,QFT,NCR(4),WS,VS,DS,GOS,WR,VR,GOR,RT
      COMMON/FIT/SIGMA,NC,NV,LAM,LEV
      DIMENSION DP(*),VC(*),DV(*),DIAG(*),PD(*),ROW(*),EVAL(*),
     &EVEC(*),AMEV(*),SAM(*),YO(*),SIGYO(*),AM(*),KI(*),IOUT(*),
     &SQSIG(2),DAMP(9)
      DIMENSION XCL(*),FCL(*),P(*),W(*),V(*),V0(*),SV(*),CS(*),
     &FOK(*),RFOK(5,*),SFOK(*),FAK(*),CAM(*),CBM(*),CCM(*),B0(*),
     &B1(*),B2(*),CP(*),ALP(2,*),BP(*),BG(*),BGP(*),CG(*),PCC(*),
     &ETHC(*),ETHCP(*),PAA(*),ETHA(*),ETHAP(*),PCCA(*),PCAA(*),
     &DCA(*),PNC(*),PNA(*),PNN(*),ZA(*),ZFA(*),ZC(*),ZFC(*),G(*),
     &GO(*),GX0(*),VRT(*),PH(*),PHV(*),XF(*),GX(*),AC(*),X(10,*),
     &XC(*),CAP(*),COP(*),XS(11,*),XCS(*),XCOS(*),XCAS(*),EXC(*),
     &WA(*),ROC(8,*),RAC(8,*),EOC(8,*),CEOC(*),REOC(5,8+NSO,*),
     &SEOC(8+NSO,*),EAC(8,*),FOC(8,*),FAC(8,*),EP(*),
     &PT(8+NSO+NSA,*),FX(*),PB(*),RZIR(*),SCPSI(*)
      DIMENSION KX(*),NF(*),ICA(*),IAN(*),INU(*),
     &NZA(*),NZPO(*),NIDT(*),NCOP(8,*),NCAP(7,*),NPSI(*)
      CHARACTER*60 LTITLE
CC
      DATA ONE,ZERO/1.0D+0,0.0D+0/
      DATA DAMP/1.0D0,0.63D0,0.40D0,0.25D0,0.16D0
     &  ,0.10D0,0.063D0,0.040D0,0.025D0/
CC
      WRITE (12,500) LTITLE
      WRITE (9,500) LTITLE
      NXP=8 
      IPCA=0
      WRITE (12,525)
      WRITE (9,525)
      WRITE (12,575) NOF
      WRITE (9,575) NOF
      IF (NC.EQ.0) GOTO 150
CC    COUNT VARIABLES
CC
130   N=0
      DO 140 I=1,NP
      IF (KI(I)) 135,140,135
135   N=N+1
      IOUT(N)=I
140   CONTINUE
      NV=N
      WRITE (12,515) NV
      WRITE (9,515) NV
      IF(NV.EQ.0) WRITE(9,485)
485   FORMAT(1X,'THERE ARE NO PARAMETERS TO BE REFINED'/)
      IF(NV.EQ.0) THEN
        PAUSE
        STOP
      END IF
      GOTO 170
150   DO 155 I=1,NP
      KI(I)=0
      DP(I)=0
155   CONTINUE
      NV=0
CC    INITIALIZE PROBLEM
CC
CC
170   NM=(NV*(NV+1))/2
      SQSIG(1)=ZERO
CC    PUT OUT TRIAL PARAMETERS, KEY INTEGERS, AND
CC
CC    PARAMETER INCREMENTS
CC
      IF(NC.EQ.0.OR.NV.EQ.0) GOTO 177 
      WRITE (12,675)
      WRITE (9,675)
      DO 175 I=1,NP
      IF(KI(I).NE.0) WRITE (12,680) I,P(I),KI(I),DP(I)
      IF(KI(I).NE.0) WRITE (9,680) I,P(I),KI(I),DP(I)
175   CONTINUE
      WRITE(12,710)
CC    START LOOP TO PERFORM NC CYCLES AND ONE FINAL
CC
CC    CALCULATION OF Y
CC
177   ICONV=0
      NCY=NC+1
      DO 440 IC=1,NCY
CC    CLEAR ARRAYS AM AND VC EXCEPT ON LAST CYCLE
CC
      IF (IC-NCY) 180,195,195
180   DO 185 I=1,NM
185   AM(I)=ZERO
      DO 190 I=1,NV
190   VC(I)=ZERO
CC    INITIALIZE FOR CYCLE IC AND PUT OUT CAPTION
CC
CC    FOR LIST OF Y(CALC)
CC
195   SQSIG(2)=SQSIG(1)
      SIG=ZERO
      IF(NCY.GT.IC) GO TO 198
      ICP=IC-1
      J=1
      DO 193 I=1,NP
      IF(KI(I).EQ.0) GO TO 193
      POLD=P(I)-PD(J)
      SIGP=DSQRT(DIAG(J))*SQSIG(1)
      J=J+1
193   CONTINUE
       WRITE(9,595)
198   CONTINUE
CC    START LOOP THROUGH NO OBSERVATIONS
CC
      DO 280 I=1,NO
CC    ENTER USERS SUBROUTINE TO COMPUTE Y(CALC)
CC
CC    AND DERIVATIVES
CC
      IDC=0
      CALL CALC(P,W,V,V0,SV,CS,G,
     &GO,GX0,VRT,PH,PHV,XF,GX,AC,X,XS,XCS,XCOS,XCAS,XC,      
     &ROC,RAC,EOC,CEOC,REOC,SEOC,EAC,FOC,FAC,FOK,RFOK,SFOK,FAK,XCL,
     &FCL,EXC,WA,YC,B0,B1,B2,CP,ALP,BP,BG,BGP,CG,PCC,
     &ETHC,ETHCP,PAA,ETHA,ETHAP,PCCA,PCAA,DCA,PNC,PNA,PNN,CAM,CBM,CCM,                                                                        
     &ZA,ZFA,ZC,ZFC,CAP,COP,EP,PT,FX,PB,RZIR,SCPSI,KX,NZA,
     &NZPO,NCOP,NCAP,NIDT,NF,ICA,IAN,INU,I,IC,IDC,ICONV,IERROR,NPSI)
CC    OBTAIN WEIGHT AND CALCULATE QUANTITIES
CC
CC    FROM Y(OBS)-Y(CALC)
CC
      SQRTW=ONE/SIGYO(I)
      DY=YO(I)-YC
      WDY=SQRTW*DY
      SIG=SIG+WDY*WDY
CC    PUT OUT Y(CALC) AND OTHER INFORMATION
CC
CC    FOR ONE OBSERVATION
CC
      IF(IC.LT.NCY) GO TO 218
      WRITE (9,625) I,YO(I),YC,DY,SIGYO(I),WDY,(X(J,I),J=1,NXP)
      WRITE (13,625) I,YO(I),YC,DY,SIGYO(I),WDY,(X(J,I),J=1,NXP)
218   CONTINUE
CC    BY-PASS DERIVATIVE AND MATRIX SET-UP ON
CC
CC    FINAL CALC OF Y
CC
      IF (IC-NCY) 220,280,280
CC    START LOOP TO STORE AN ARRAY OF NV
CC
CC    DERIVATIVES
CC
220   J=1
      DO 250 K=1,NP
      IF (KI(K)) 250,250,235
CC    OBTAIN DERIVATIVE NUMERICALLY
CC
235   DPK=DP(K)
      IF (DPK.NE.0) GOTO 240
      WRITE(9,490) K
490   FORMAT(1X,'PARAMETER ',I2,' INCREMENT IS ZERO')
      PAUSE
      STOP
240   PSAVE=P(K)
      P(K)=PSAVE+DPK
      CALL CALC(P,W,V,V0,SV,CS,G,
     &GO,GX0,VRT,PH,PHV,XF,GX,AC,X,XS,XCS,XCOS,XCAS,XC,
     &ROC,RAC,EOC,CEOC,REOC,SEOC,EAC,FOC,FAC,FOK,RFOK,SFOK,FAK,XCL,
     &FCL,EXC,WA,YD,B0,B1,B2,CP,ALP,BP,BG,BGP,CG,PCC,
     &ETHC,ETHCP,PAA,ETHA,ETHAP,PCCA,PCAA,DCA,PNC,PNA,PNN,CAM,CBM,CCM,
     &ZA,ZFA,ZC,ZFC,CAP,COP,EP,PT,FX,PB,RZIR,SCPSI,KX,NZA,
     &NZPO,NCOP,NCAP,NIDT,NF,ICA,IAN,INU,I,IC,J,ICONV,IERROR,NPSI)
      DV(J)=SQRTW*(YD-YC)/DPK
      P(K)=PSAVE
245   J=J+1
250   CONTINUE
CC    END LOOP TO OBTAIN DERIVATIVES
CC
CC    START LOOP TO STORE MATRIX AND VECTOR.
CC
C     1604 OR GLS STORAGE SCHEME IS REVERSE
C
C     OF 7090 OR GLS
C
      JK=1
      DO 275 J=1,NV
CC
255   TEMP=DV(J)
      IF (TEMP) 265,260,265
CC    BY-PASS IF DERIVATIVE IS ZERO
CC
260   JK=JK+NV+1-J
      GO TO 275
265   DO 270 K=J,NV
      AM(JK)=AM(JK)+TEMP*DV(K)
      JK=JK+1
270   CONTINUE
      VC(J)=VC(J)+TEMP*WDY
275   CONTINUE
CC    END LOOP TO STORE MATRIX AND
CC
CC    VECTOR
CC
280   CONTINUE
CC    END LOOP THROUGH NO OBSERVATIONS
CC
CC    COMPUTE AND PUT OUT AGREEMENT FACTORS
CC
      SQSIG(1)=DSQRT(SIG/DBLE(FLOAT(MAX0(NOF-NV,1))))
      SIGSV=SIG
      SQSGSV=SQSIG(1)
      SIGMA=SQSIG(1)
      WRITE (12,630) IC,SIG,SQSIG(1)
      WRITE (9,630) IC,SIG,SQSIG(1)
CC    BY-PASS MATRIX INVERSION AND PARAMETER OUTPUT
CC
CC    ON FINAL CYCLE
CC
      IF (IC-NCY) 285,445,445
CC    SAVE MATRIX FOR POSSIBLE PRINCIPAL COMPONENT
CC
CC    ANALYSIS
CC
285   JK=0
      DO 295 J=1,NV
      DO 290 K=J,NV
      JK=JK+1
      JKU=J+(K-1)*NV
      JKL=K+(J-1)*NV
      EVEC(JKU)=AM(JK)
      EVEC(JKL)=AM(JK)
290   CONTINUE
295   CONTINUE
CC    START LOOP TO TEST FOR ZERO DIAGONAL ELEMENT
CC
      ISING=0
      II=1
      IID=NV
      DO 310 I=1,NV
      IF (AM(II)) 305,300,305
300   ISING=1
      WRITE (12,640) IOUT(I)
      WRITE (9,640) IOUT(I)
305   II=II+IID
      IID=IID-1
310   CONTINUE
CC    END LOOP TO TEST FOR ZERO DIAGONAL
CC
CC    ELEMENT
CC
CC    PRINCIPAL COMPONENT ANALYSIS IF ZERO DIAGONAL
CC
CC    FOUND
CC
      IF (ISING) 330,315,330
CC    ENTER SUBROUTINE TO REPLACE MATRIX
CC
CC    WITH INVERSE
CC
315   CALL MATINV(AM,NV,ISING)
      IF (ISING) 325,320,325
320   IF (IPCA) 330,365,330
CC    SINGULAR MATRIX FOUND
CC
325   WRITE (12,650)
      WRITE (9,650)
CC    PRINT MATRIX
CC
330   WRITE (12,600) IC,(IOUT(I),I=1,NV)
      WRITE (9,600) IC,(IOUT(I),I=1,NV)
      K2=0
      DO 335 I=1,NV
      K1=K2+1
      K2=K2+NV
      WRITE (12,605) IOUT(I),(EVEC(K),K=K1,K2)
      WRITE (9,605) IOUT(I),(EVEC(K),K=K1,K2)
335   CONTINUE
CC    COMPUTE AND PRINT SCALE FACTORS
CC
      CALL MSTR(EVEC,SAM,NV,0,2)
      DO 345 I=1,NV
      IF (SAM(I)) 345,345,340
340   SAM(I)=ONE/DSQRT(SAM(I))
345   CONTINUE
      WRITE (12,620) (IOUT(I),I=1,NV)
      WRITE (9,620) (IOUT(I),I=1,NV)
      WRITE (12,645) (SAM(I),I=1,NV)
      WRITE (9,645) (SAM(I),I=1,NV)
CC    SCALE MATRIX
CC
      IJ=0
      DO 355 J=1,NV
      DO 350 I=1,NV
      IJ=IJ+1
      EVEC(IJ)=EVEC(IJ)*SAM(I)*SAM(J)
350   CONTINUE
355   CONTINUE
CC    PRINCIPAL COMPONENT ANALYSIS
CC
      CALL MSTR(EVEC,AMEV,NV,0,1)
      CALL EIGEN(AMEV,EVEC,NV,0)
      CALL MSTR(AMEV,EVAL,NV,1,2)
CC    PRINT EIGENVALUES AND EIGENVECTORS
CC
      WRITE (12,610) IC,(IOUT(I),I=1,NV)
      WRITE (9,610) IC,(IOUT(I),I=1,NV)
      K2=0
      DO 360 I=1,NV
      K1=K2+1
      K2=K2+NV
      WRITE (12,615) I,EVAL(I),(EVEC(K),K=K1,K2)
      WRITE (9,615) I,EVAL(I),(EVEC(K),K=K1,K2)
360   CONTINUE
CC    TERMINATE IF MATRIX IS SINGULAR
CC
      IF (ISING.NE.0) WRITE(12,1111)
1111  FORMAT(1X,'@')
      IF (ISING.NE.0) RETURN
CC    START LOOP FOR MATRIX VECTOR MULTIPLICATION
CC    FOR
CC    PARAMETER
CC
CC    CHANGES
CC
365   DO 390 I=1,NV
      PDI=ZERO
      IJ=I
      IJD=NV-1
      DO 385 J=1,NV
      PDI=PDI+AM(IJ)*VC(J)
      IF (J-I) 370,375,380
370   IJ=IJ+IJD
      IJD=IJD-1
      GO TO 385
CC    SAVE DIAGONAL ELEMENTS
CC
CC    OF INVERSE MATRIX
CC
375   DIAG(I)=AM(IJ)
380   IJ=IJ+1
385   CONTINUE
      PD(I)=PDI
      SIG=SIG-PDI*VC(I)
390   CONTINUE
CC    END LOOP FOR MATRIX VECTOR MULTIPLICATION
CC
CC    RECOMPUTE AGREEMENT FACTOR USING
CC
CC    MODIFIED SIG
CC
      SQSIG(1)=DSQRT(DMAX1(SIG/ FLOAT(MAX0(NOF-NV,1)),ZERO))
CC    PUT OUT CAPTION FOR LIST OF CORRECTED
CC
CC    PARAMETERS
CC
      WRITE (12,655) IC
      WRITE (9,655) IC
CC    START LOOP TO CORRECT AND PUT OUT
CC
CC    PARAMETERS
CC
      J=1
      DO 405 I=1,NP
      IF(KI(I).EQ.0) GO TO 405
      POLD=P(I)
      KII=KI(I)
      PD(J)=PD(J)*DAMP(KII)
      P(I)=POLD+PD(J)
      SIGP=DSQRT(DIAG(J))*SQSIG(1)
      WRITE (12,665) I,POLD,PD(J),P(I),SIGP
      WRITE (9,665) I,POLD,PD(J),P(I),SIGP
      J=J+1
405   CONTINUE
CC    END LOOP TO CORRECT AND PUT OUT
CC
CC    PARAMETERS
CC
CC    PUT OUT ESTIMATED AGREEMENT FACTORS
CC
      WRITE (9,635) IC,SIG,SQSIG(1)
      WRITE (12,635) IC,SIG,SQSIG(1)
CC    ENTER USERS SUBROUTINE TO TEST AND
CC
CC    MODIFY PARAMETERS
CC
CC    OR END JOB
CC
      ISTOP='CONTINUE'
      CALL TEST(P,PD,KI,ICONV,ISTOP)
CC
CC    TERMINATE REFINEMENT IF INDICATED BY USERS
CC
CC    SUBROUTINE TEST
CC
      IF (ISTOP.EQ.'CONTINUE') GO TO 440
      WRITE (12,670) ISTOP
      WRITE (9,670) ISTOP
      NCY=IC+1
      NC=IC
440   CONTINUE
CC    END LOOP THROUGH NC CYCLES AND
CC
CC    FINAL CALC OF Y
CC
CC    TERMINATE JOB
CC
445   IF (NC.LE.0) WRITE(12,1111)
      IF (NC.LE.0) RETURN
CC    CALCULATE AND PUT OUT CORRELATION MATRIX
CC
450   WRITE (9,700)
      WRITE (12,700)
      DO 455 I=1,NV
      DIAG(I)=ONE/DSQRT(DIAG(I))
455   CONTINUE
      IJ=1
      DO 470 I=1,NV
      DO 460 J=1,NV
      ROW(J)=ZERO
460   CONTINUE
      DO 465 J=I,NV
      ROW(J)=AM(IJ)*DIAG(I)*DIAG(J)
      IJ=IJ+1
465   CONTINUE
      WRITE (12,705) IOUT(I),(ROW(J),J=1,NV)
      WRITE (9,705) IOUT(I),(ROW(J),J=1,NV)
470   CONTINUE
      SIGMA=SQSIG(1)
      WRITE(12,1111)
      DO 475 I=1,NP
      IF(KI(I).NE.0) WRITE(18,1112) I,P(I)
1112  FORMAT(1X,I3,1PE12.4)
475   CONTINUE
      RETURN
CC    FORMAT STATEMENTS
500   FORMAT(//1X,A60/)
515   FORMAT(1X,'NUMBER OF PARAMETERS TO BE VARIED IS',I3)
525   FORMAT(1X,'PRINCIPAL COMPONENT ANALYSIS IN CASE OF ERROR'/)
575   FORMAT(1X,'NUMBER OF OBSERVATIONS READ IS',I4/)
590   FORMAT(1X,'CALCULATED Y BASED ON PARAMETERS BEFORE CYCLE',I2)
595   FORMAT(1X,' PT     Y(OBS)    Y(CALC)   OBS-CALC     SIG(O)',
     &'  (O-C)/SIG       X(I)'/1H )
600   FORMAT(1X,'MATRIX FOR CYCLE',I3/'0',3X,10I11/(4X,10I11))
605   FORMAT(1X,I3,10E11.3/(4X,10E11.3))
610   FORMAT(1X,'PRINCIPAL COMPONENTS FOR CYCLE',I3/1X,'  I  EIGVAL(I)',
     &5X,10I10/(20X,10I10))
615   FORMAT(1X,I3,E11.3,5X,10F10.6/(20X,10F10.6))
620   FORMAT(1X,'SCALE FACTORS APPLIED TO MATRIX'/'0',10I11/(1X,10I11))
625   FORMAT(1X,I3,1P7E11.3/(14X,6E11.3))
630   FORMAT(/1X,'AGREEMENT FACTORS BASED ON PARAMETERS BEFORE CYCLE',
     &I2/1X,'SUM(W*(O-C)**2) IS ',1PE11.3/1X,
     &'SQRTF(SUM(W*(O-C)**2)/(NO-NV)) IS ',E12.4)
635   FORMAT(/1X,'ESTIMATED AGREEMENT FACTORS BASED ON PARAMETERS',
     &' CYCLE',I2/1X,'SUM(W*(O-C)**2) IS ',1PE11.3/
     &1X,'SQRTF(SUM(W*(O-C)**2)/(NO-NV)) IS ',E12.4)
640   FORMAT(1X,' MATRIX HAS A ZERO DIAGONAL ELEMENT CORRESPONDING TO PAR
     &AMETER',I3)
645   FORMAT(1X,10E11.3/(1X,10E11.3))
650   FORMAT(1X,' SINGULARITY RETURN FROM MATRIX INVERTER')
655   FORMAT(/1X,'PARAMETERS AFTER LEAST SQUARES CYCLE',I2,//1X,
     &'           OLD         CHANGE        NEW          ERROR',/)
665   FORMAT(' ',I3,1P4E14.5)
670   FORMAT(1X,'SUBROUTINE TEST INDICATES THAT JOB IS TO BE TERMINATED
     & FOR REASON: ',A10)
675   FORMAT(1X,'STARTING VALUES'/
     &/1X,'  I      P(I)      KI(I)      DP(I)'/1X,' ')
680   FORMAT(1X,I3,3X,1PE12.5,3X,I4,3X,E11.3)
700   FORMAT(/1X,'CORRELATION MATRIX')
705   FORMAT(1X,I3,10F9.4/(1X,3X,10F9.4))
710   FORMAT(/)
      END
C
C
C
      SUBROUTINE EIGEN(A,R,N,MV)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(*),R(*)
      RETURN
      END
C
C
      SUBROUTINE MATINV(AM,N,NFAIL)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      DIMENSION AM(*)
C     ********** SEGMENT 1 OF CHOLESKI INVERSION **********
C
CC    ***** FACTOR MATRIX INTO LOWER TRIANGLE X TRANSPOSE
CC
CC    *****
CC
      DATA ONE,ZERO/1.0D+0,0.0D+0/
      K=1
      IF (N-1) 2,5,10
2     NFAIL=K
      GOTO 105
5     AM(1)=ONE/AM(1)
      GO TO 100
CC    ***** LOOP M OF A(L,M) *****
CC
10    DO 55 M=1,N
      IMAX=M-1
CC    ***** LOOP L OF A(L,M) *****
CC
      DO 50 L=M,N
      SUMA=ZERO
      KLI=L
      KMI=M
      IF (IMAX) 25,25,15
C     *****SUM OVER I=1,M-1 A(L,I)*A(M,I) *****
C
15    DO 20 I=1,IMAX
      SUMA=SUMA+AM(KLI)*AM(KMI)
      J=N-I
      KLI=KLI+J
20    KMI=KMI+J
CC    *****TERM=C(L,M)-SUM *****
CC
25    TERM=AM(K)-SUMA
      IF (L-M) 30,30,45
30    IF (TERM) 40,40,35
CC    ***** A(M,M)=SQRT(TERM) *****
CC
35    DENOM=DSQRT(TERM)
      AM(K)=DENOM
      GO TO 50
40    NFAIL= K
      GO TO 105
CC    ***** A(L,M)=TERM/A(M,M) *****
CC
45    AM(K)=TERM/DENOM
50    K=K+1
55    CONTINUE
C     ********** SEGMENT 2 OF CHOLESKI INVERSION **********
C
CC    *****INVERSION OF TRIANGULAR MATRIX*****
CC
60    AM(1)=ONE/AM(1)
      KDM=1
CC    ***** STEP L OF B(L,M) *****
CC
      DO 75 L=2,N
      KDM=KDM+N-L+2
CC    ***** RECIPROCAL OF DIAGONAL TERM *****
CC
      TERM = ONE/AM(KDM)
      AM(KDM)=TERM
      KMI=0
      KLI=L
      IMAX=L-1
CC    ***** STEP M OF B(L,M) *****
CC
      DO 70 M=1,IMAX
      K=KLI
CC    ***** SUM TERMS *****
CC
      SUMA=ZERO
      DO 65 I=M,IMAX
      II=KMI+I
      SUMA=SUMA-AM(KLI)*AM(II)
65    KLI=KLI+N-I
CC    ***** MULT SUM * RECIP OF DIAGONAL *****
CC
      AM(K)=SUMA*TERM
      J=N-M
      KLI=K+J
70    KMI=KMI+J
75    CONTINUE
C     ********** SEGMENT 3 OF CHOLESKI INVERSION **********
C
CC    *****PREMULTIPLY LOWER TRIANGLE BY TRANSPOSE*****
CC
80    K=1
      DO 95 M=1,N
      KLI=K
      DO 90 L=M,N
      KMI=K
      IMAX=N-L+1
      SUMA=ZERO
      DO 85 I=1,IMAX
      SUMA=SUMA+AM(KLI)*AM(KMI)
      KLI=KLI+1
85    KMI=KMI+1
      AM(K)=SUMA
90    K=K+1
95    CONTINUE
100   NFAIL=0
105   RETURN
      END
C
C
C
      SUBROUTINE MSTR(A,R,N,MSA,MSR)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
C     THIS IS ORNL F01042 OF 1167
C
CC    ..................................................................
CC
CC    SUBROUTINE MSTR
CC
CC    PURPOSE
CC
CC    CHANGE STORAGE MODE OF A MATRIX
CC
CC    USAGE
CC
CC    CALL MSTR(A,R,N,MSA,MSR)
CC
CC    DESCRIPTION OF PARAMETERS
CC
CC    A - NAME OF INPUT MATRIX
CC
CC    R - NAME OF OUTPUT MATRIX
CC
CC    N - NUMBER OF ROWS AND COLUMNS IN A AND
CC
CC    R
CC
CC    MSA - ONE DIGIT NUMBER FOR STORAGE MODE
CC
CC    OF MATRIX A
CC
CC    0 - GENERAL
CC
C     1 - SYMMETRIC
C     2 - DIAGONAL
C
CC    MSR - SAME AS MSA EXCEPT FOR MATRIX R
CC
CC    REMARKS
CC
CC    MATRIX R CANNOT BE IN THE SAME LOCATION
CC
CC    AS MATRIX A
CC
CC    MATRIX A MUST BE A SQUARE MATRIX
CC
CC    SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
CC
CC    LOK (FORMERLY LOC)
CC
CC    METHOD
CC
CC    MATRIX A IS RESTRUCTURED TO FORM MATRIX
CC
CC    R.
CC
CC    MSA MSR
CC
CC    0 0 MATRIX A IS MOVED TO MATRIX R
CC
C     0 1 THE UPPER TRIANGLE ELEMENTS OF
C
C     A GENERAL MATRIX
C
CC    ARE USED TO FORM A SYMMETRIC
CC
CC    MATRIX
CC
C     0 2 THE DIAGONAL ELEMENTS OF A GENERAL
C
C     MATRIX ARE USED
C
CC    TO FORM A DIAGONAL MATRIX
CC
C     1 0 A SYMMETRIC MATRIX IS EXPANDED
C
C     TO FORM A GENERAL
C
CC    MATRIX
CC
C     1 1 MATRIX A IS MOVED TO MATRIX R
C
C     1 2 THE DIAGONAL ELEMENTS OF A SYMMETRIC
C     MATRIX ARE
C
CC    USED TO FORM A DIAGONAL MATRIX
CC
C     2 0 A DIAGONAL MATRIX IS EXPANDED
C
C     BY INSERTING MISSING
C
CC    ZERO ELEMENTS TO FORM A GENERAL
CC
CC    MATRIX
CC
C     2 1 A DIAGONAL MATRIX IS EXPANDED
C
C     BY INSERTING MISSING
C
CC    ZERO ELEMENTS TO FORM A SYMMETRIC
CC
CC    MATRIX
CC
C     2 2 MATRIX A IS MOVED TO MATRIX R
C
CC
CC    ..................................................................
CC
      DIMENSION A(*),R(*)
      DATA ONE,ZERO/1.0D+0,0.0D+0/
      DO 25 I=1,N
      DO 25 J=1,N
CC    IF R IS GENERAL, FORM ELEMENT
CC
      IF (MSR) 5,10,5
CC    IF IN LOWER TRIANGLE OF SYMMETRIC OR DIAGONAL
CC
CC    R, BYPASS
CC
5     IF (I-J) 10,10,25
10    CALL LOK(I,J,IR,N,N,MSR)
CC    IF IN UPPER AND OFF DIAGONAL OF DIAGONAL
CC
CC    R, BYPASS
CC
      IF (IR) 25,25,15
CC    OTHERWISE, FORM R(I,J)
CC
15    R(IR)=ZERO
      CALL LOK(I,J,IA,N,N,MSA)
CC    IF THERE IS NO A(I,J), LEAVE R(I,J) AT ZERO
CC
      IF (IA) 25,25,20
20    R(IR)=A(IA)
25    CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE LOK(I,J,IR,N,M,MS)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
CC    ..................................................................
CC
CC    SUBROUTINE LOK (FORMERLY LOC)
CC
CC    PURPOSE
CC
CC    COMPUTE A VECTOR SUBSCRIPT FOR AN ELEMENT
CC
CC    IN A MATRIX OF
CC
CC    SPECIFIED STORAGE MODE
CC
CC    USAGE
CC
CC    CALL LOK (I,J,IR,N,M,MS)
CC
CC    DESCRIPTION OF PARAMETERS
CC
CC    I - ROW NUMBER OF ELEMENT
CC
CC    J - COLUMN NUMBER OF ELEMENT
CC
CC    IR - RESULTANT VECTOR SUBSCRIPT
CC
CC    N - NUMBER OF ROWS IN MATRIX
CC
CC    M - NUMBER OF COLUMNS IN MATRIX
CC
CC    MS - ONE DIGIT NUMBER FOR STORAGE MODE
CC
CC    OF MATRIX
CC
CC    0 - GENERAL
CC
C     1 - SYMMETRIC
C     2 - DIAGONAL
C
CC
CC    REMARKS
CC
CC    NONE
CC
CC    SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
CC
CC    NONE
CC
CC    METHOD
CC
CC    MS=0 SUBSCRIPT IS COMPUTED FOR A MATRIX
CC
CC    WITH N*M ELEMENTS
CC
CC    IN STORAGE (GENERAL MATRIX)
CC
C     MS=1 SUBSCRIPT IS COMPUTED FOR A MATRIX
C
C     WITH N*(N+1)/2 IN
C
CC    STORAGE (UPPER TRIANGLE OF SYMMETRIC
CC
CC    MATRIX). IF
CC
CC    ELEMENT IS IN LOWER TRIANGULAR
CC
CC    PORTION, SUBSCRIPT IS
CC
CC    CORRESPONDING ELEMENT IN UPPER
CC
CC    TRIANGLE.
CC
C     MS=2 SUBSCRIPT IS COMPUTED FOR A MATRIX
C
C     WITH N ELEMENTS
C
CC    IN STORAGE (DIAGONAL ELEMENTS OF
CC
CC    DIAGONAL MATRIX).
CC
CC    IF ELEMENT IS NOT ON DIAGONAL (AND
CC
CC    THEREFORE NOT IN
CC
CC    STORAGE), IR IS SET TO ZERO.
CC
CC    ..................................................................
CC
      IX=I
      JX=J
      IF (MS-1) 5,10,25
5     IRX=N*(JX-1)+IX
      GO TO 35
10    IF (IX-JX) 15,20,20
15    IRX=IX+(JX*JX-JX)/2
      GO TO 35
20    IRX=JX+(IX*IX-IX)/2
      GO TO 35
25    IRX=0
      IF (IX-JX) 35,30,35
30    IRX=IX
35    IR=IRX
      RETURN
      END
C
C
C
      SUBROUTINE TEST(P,PD,KI,ICONV,ISTOP)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      CHARACTER*10 ISTOP
      COMMON/LGP/NP
      COMMON/MD/NSO,NCBO,IREF,NPH,NPA,NSA,NIO,NIA,NCA,NAN,NNU,NUP,NUM
      COMMON/FIT/SIGMA,NC,NV,LAM,LEV     
      DIMENSION P(*),PD(*),KI(*)
      DATA ONE,ZERO/1.0D+0,0.0D+0/
      PDSV=ZERO
      NAP=1+NPH+NPA
      J=0
      DO 10 I=1,NP
      IF(KI(I).EQ.0) GO TO 10
      J=J+1
      IF(NSO.EQ.0) GOTO 14
      DO 12 K=1,NSO
      L=NAP*(K-1)+1
      IF(I.EQ.L) GOTO 22
12    CONTINUE
14    IF(NSA.EQ.0) GOTO 18
      DO 16 K=1,NSA
      L=NAP*NSO+2*(NPH+NPA)
      IF(I.GT.L.AND.I.LE.(L+NSA)) GOTO 22
16    CONTINUE
C  PARAMETER IS NOT LOGARITHMIC
18    IF(DABS(P(I)-PD(J)).LT.1) DPF=PD(J)
      IF(DABS(P(I)-PD(J)).GE.1) DPF=PD(J)/(P(I)-PD(J))
      ABSPD=DABS(DPF)
      IF(ABSPD.GT.PDSV) PDSV=ABSPD
      GOTO 10
C  PARAMETER IS LOGARITHMIC
22    DPF=PD(J)
      ABSPD=DABS(DPF)
      IF(ABSPD.GT.PDSV) PDSV=ABSPD
10    CONTINUE
      IF(PDSV.GT.0.1D+0) ICONV=0
      IF(PDSV.LE.0.1D+0) ICONV=1
      IF(PDSV.LE.ONE) GO TO 40
      WRITE(9,100) PDSV
      WRITE(12,100) PDSV
100   FORMAT(//1X,'REDUCING PARAMETER CHANGES BY A FACTOR OF '1PE9.2/)
      J=0
      DO 50 I=1,NP
      IF(KI(I).EQ.0) GO TO 50
      J=J+1
      P(I)=P(I)-PD(J)+PD(J)/PDSV
50    CONTINUE
      RETURN
40    IF(PDSV.LE.1.0D-7) ISTOP='CONVERGED'
C
C      WRITE(9,1000) ICONV
C      WRITE(12,1000) ICONV
C1000  FORMAT(1X,'ICONV =',I3)
C
      RETURN
      END
C
C
C
C  SUBROUTINE THAT CONTROLS THE CALCULATION OF THE OBSERVED QUANTITY
C
C  
C                              GLOSSARY
C
C     CS(I) = EQUILIBRIUM CONCENTRATION OF SPECIES I BEING CALCD.
C       DEN = DENSITY OF THE AQUEOUS PHASE (G/CM3). 
C      DEXC = INITIAL VALUE OF INCREMENTS USED IN SIMPLEX CALCN. 
C        EP = ARRAY USED BY SMPLX
C       EXC = ARRAY USED BY SMPLX
C    FCL(I) = RESIDUAL IN EQUATION I, TO BE TAKEN TO ZERO. 
C      FTOL = ERROR LIMIT USED IN HYBRD2.
C        FX = ARRAY USED BY SMPLX
C      G(I) = ACTIVITY COEFFICIENT OF SPECIES I
C       IAO = INDICATES PRESENCE OF (0) NONAQUEOUS PHASE ONLY, (1) BOTH
C              AQUEOUS AND NONAQUEOUS PHASES, OR (2) AQ. PHASE ONLY.
C        IC = REFINEMENT CYCLE NUMBER.
C       ICN = INDICATOR OF CONC. SCALE; MOLARITY (1), MOLALITY (2),
C              OR MOLES/KG SOLN.
C       IDC = IF > 0, INDICATES CALCULATION OF A PARTIAL DERIVATIVE.
C      IDEN = INDICATOR OF AQ.PROD.SPECIES FORMATION IN THE ABSENCE 
C              OF HA OR B.
C       IEH = INDICATOR OF GIVEN INITIAL (0) OR EQUIL.AQ.ACIDITY (1),
C              OR GIVEN EQUIL. AQ. H+ CONC. (>1). 
C       IEQ = INDICATOR OF GIVEN INITIAL (0) OR EQUIL. (1) AQ.COMPN.
C       IER = ERROR CONDITION ON RETURN FROM HYBRD2.
C    IERROR = ERROR CONDITION ON RETURN TO ORGLS.
C       IPT = THE NO. OF THE CURRENT DATA POINT.
C        IY = INDICATOR OF KIND OF OBSERVATION TO BE CALCULATED..
C        KX = ARRAY USED BY SMPLX
C       LWA = LENGTH OF ARRAY WA 
C    MAXFEV = NUMBER OF CALLS TO ENTFCN FROM HYBRD2 (MAX.ON ENTRY
C              AND ACTUAL ON RETURN FROM HYBRD2).
C      NASP = NUMBER OF POSSIBLE AQ.PROD.SPECIES NOT INVOLVING HA AND/OR B
C      NCAL = FLAG INDICATING THE KIND OF EQUILIBRIUM CALCULATION:
C              0, NONE NEEDED; 1, BOTH AQ.AND NONAQ.PHASES PRESENT;
C              2, ONLY NONAQ.PHASE PRESENT; 3, ONLY AQ.PHASE PRESENT
C NCAP(I,J) = STOICHIOMETRY COEFF.OF REACTANT I IN AQ. PROD.SPEC.J. 
C    NCR(I) = THE DEGREE OF ASSOCIATION OF HA (1) OR B (2) 
C      NCSM = MAXIMUM NUMBER OF CYCLES USED IN SIMPLEX CALCULATION. 
C       NEQ = THE NUMBER OF SPEC.CONCNS.TO BE CALCULATED 
C        NF = A LABEL DESIGNATING EACH CONCENTRATION TO BE DETMD.
C       NIA = INDICATOR OF ASSUMED UNIT VALUES FOR AQ.ACT.COEFS.
C   NIDT(I) = SIX DIGIT INTEGER SPECIFYING DATA TYPE
C        NO = NUMBER OF DATA POINTS
C       NSA = NUMBER OF POSSIBLE AQ.PROD.SPECIES 
C       NSO = NUMBER OF POSSIBLE NONAQ.PROD.SPECIES 
C       NSX = NUMBER OF CYCLES USED IN SIMPLEX CALCULATION.
C        PB = ARRAY USED BY SMPLX
C        PT = ARRAY USED BY SMPLX
C      RMBA = MATERIAL BALANCE CORRECTION FACTOR FOR THE AQUEOUS PHASE
C      RMBO = MATERIAL BALANCE CORRECTION FACTOR FOR THE NONAQUEOUS PHASE
C      SMIN = VALUE OF TFCL ABOVE WHICH SIMPLEX CALCN.IS REQUIRED.
C      SMTL = CONVERGENCE LIMIT USED IN SIMPLEX 
C       TFC = SUM OF SQUARES OF RESIDUALS FCL 
C      TFCL = SQUARE ROOT OF TFC  
C      TFCS = CURRENT LOWEST VALUE OF TFC 
C      TSAI = INITIAL AMOUNT OF AQUEOUS SOLUTE (CM3/L,G/KG, OR G/KG SLN.)
C      TSOI = INITIAL AMOUNT OF NONAQUEOUS SOLUTE (CM3/L,G/KG, OR G/KG SLN.)
C      TVSA = APPARENT VOLUME OF SOLUTE (CM3/L.SLN.,CM3/KG SLV.,CM3/KG SLN.)
C      V(I) = MOLAR VOLUME OF SPECIES (CM3/MOL) 
C      VWPA = AMOUNT OF WATER RELEASED BY PRODUCT SPECIES FORMATION IN THE
C              AQUEOUS PHASE (CM3/L,G/KG, OR G/KG SLN.)
C      VWPO = AMOUNT OF WATER RELEASED BY PRODUCT SPECIES FORMATION IN THE
C              NONAQUEOUS PHASE (CM3/L AQ.SLN.,G/KG H2O, OR G/KG AQ.SLN.)
C      W(I) = MOLECULAR WEIGHT OF SPECIES (G/MOL)
C        WA = SCRATCH ARRAY USED IN HYBRD2
C    X(I,J) = THE GIVEN CONC. OF REACTANT I FOR POINT J. 
C     XC(I) = THE EST.LOG.EQUILIBRIUM CONCENTRATION OF SPECIES I. 
C      XCAS = A LINEAR ARRAY WHICH STORES THE XC VALUES OF ALL SPECIES
C              CALCULATED PREVIOUSLY WHEN ONLY THE AQUEOUS PHASE IS
C              CONSIDERED.
C      XCOS = A LINEAR ARRAY WHICH STORES THE XC VALUES OF ALL SPECIES
C              CALCULATED PREVIOUSLY WHEN ONLY THE NONAQUEOUS PHASE IS
C              PRESENT.
C       XCS = A LINEAR ARRAY WHICH STORES THE XC VALUES OF
C              ALL SPECIES CALCULATED PREVIOUSLY WHEN BOTH PHASES ARE
C              PRESENT.
C      XTOL = ERROR LIMIT USED IN HYBRD2 
C
C
C            CALCULATION OF SPECIE CONCENTRATIONS AT EQUILIBRIUM
C
C  THE SPECIE CONCENTRATIONS EVALUATED FOR A DATA POINT ARE STORED IN ARRAY
C   XC IN THE FOLLOWING ORDER:
C
C              XC(1) = LOG. EQ. NONAQ. CONC. OF HA REACTANT SPECIES
C                      AS MONOMER
C              XC(2) = LOG. EQ. NONAQ. CONC. OF B REACTANT SPECIES
C                      AS MONOMER
C              XC(3) = LOG. EQ. AQ. CONC. OF M\ZM+
C              XC(4) = LOG. EQ. AQ. CONC. OF N\ZN+
C              XC(5) = LOG. EQ. AQ. CONC. OF X\ZX+
C              XC(6) = LOG. EQ. AQ. CONC. OF Y\ZY-
C              XC(7) = LOG. EQ. AQ. CONC. OF H\+
C              XC(8) = LOG. EQ. AQ. CONC. OF OH\-
C              XC(9) = LOG. EQ. CONC. OF NONAQ. PROD. SP. NO. 1
C          XC(8+NSO) = LOG. EQ. CONC. OF NONAQ. PROD. SP. NO. NSO
C          XC(9+NSO) = LOG. EQ. CONC. OF AQU. PROD. SP. NO. 1
C      XC(8+NSO+NSA) = LOG. EQ. CONC. OF AQU. PROD. SP. NO. NSA
C
C  NOT ALL THESE QUANTITIES NEED BE CALCULATED FOR ALL POINTS:  FOR SOME
C   DATA POINTS, SOME REACTANTS MAY NOT BE PRESENT; SOME REACTANTS MAY NOT
C   BE NEEDED TO PRODUCE THE PRODUCT SPECIES ASSUMED; ONE PHASE MAY BE 
C   ABSENT AND HENCE SOME REACTANTS WILL NOT BE PRESENT AND SOME PRODUCT
C   SPECIES MAY NOT FORM.
C
C  FOR A GIVEN CASE, THE NEQ QUANTITIES IN THIS LIST TO BE CALCULATED
C   ARE ENTERED INTO AN ARRAY XCL SEQUENTIALLY, ASCENDING IN THE ORDER
C   SHOWN ABOVE.  THEIR POSITION IN THIS ARRAY IS STORED IN AN INTEGER
C   ARRAY NF.  THE ARRAY FCL CONTAINS THE RESIDUALS OF THE NEQ EQUATIONS
C   TO BE SATISFIED IN SOLVING FOR THE EQUILIBRIUM CONCENTRATIONS IN XCL.
C
C  THERE ARE THREE KINDS OF EQUILIBRIUM CALCULATIONS, DEFINED AS FOLLOWS:
C
C    NCAL=1 (FLAG IAOC=1)
C
C        A CALCULATION EMPLOYED WHEN BOTH THE NONAQUEOUS PHASE AND THE AQUEOUS
C        PHASES ARE PRESENT. THE RESULTING SPECIE CONCENTRATIONS FOR EACH
C        POINT ARE STORED IN THE ARRAY XCS IN THE SAME ORDER THEY ARE
C        STORED IN ARRAY XC.  THIS ARRAY IS DIMENSIONED XCS(NO*(8+NSO+NSA)).
C
C    NCAL=2 (FLAG IOC=1)
C
C        A CALCULATION EMPLOYED WHEN ONLY THE NONAQUEOUS PHASE IS PRESENT. IN
C        THIS CASE IT IS NECESSARY TO ASSUME THAT AN INFINITESIMALLY SMALL
C        AMOUNT OF AN AQUEOUS PHASE IS ALSO PRESENT BECAUSE OF THE WAY
C        FORMATION CONSTANTS OF PRODUCT SPECIES ARE DEFINED. FOR EACH POINT,
C        THE RESULTING SPECIE CONCENTRATIONS, INCLUDING THOSE IN THE TRACE
C        AMOUNT OF AQUEOUS PHASE, ARE STORED IN THE ARRAY XCOS IN THE SAME
C        ORDER THEY ARE STORED IN ARRAY XC.  THIS ARRAY IS DIMENSIONED
C        XCOS(NO*(8+NSO+NSA)).
C
C    NCAL=3 (FLAG IAC=1)
C
C        A CALCULATION EMPLOYED WHEN ONLY THE AQUEOUS PHASE IS PRESENT, OR
C        WHEN AN INITIAL CALCULATION OF EQUIL.CONDITIONS IN THE AQUEOUS
C        PHASE BEFORE CONTACT WITH THE NONAQUEOUS PHASE PHASE IS NEEDED FOR
C        CALCULATING MATERIAL BALANCES.  THE RESULTING SPECIE CONCENTRATIONS
C        ARE STORED IN THE ARRAY XCAS IN THE SAME ORDER THEY ARE STORED IN
C        ARRAY XC, EXCEPT THAT THE REACTANTS HA AND B, AND ANY PRODUCT
C        SPECIES IN THE NONAQUEOUS PHASE ARE EXCLUDED.  THIS ARRAY IS THERE-
C        FORE DIMENSIONED XCAS(NO*(6+NSA)).
C
C
C  THE QUANTITIES THAT DETERMINE WHICH CALCULATIONS MUST BE PERFORMED FOR
C   A DATA POINT ARE:
C
C       IAO = INDICATES PRESENCE OF (0) NONAQUEOUS PHASE ONLY, (1) BOTH
C              AQUEOUS AND NONAQUEOUS PHASES, OR (2) AQ. PHASE ONLY.
C      NASP = NUMBER OF POSSIBLE AQ.PROD.SPECIES NOT INVOLVING HA AND/OR B
C       IEQ = INDICATOR OF GIVEN INITIAL (0) OR EQUIL. (1) AQ.COMPN.
C       ICN = INDICATOR OF CONC. SCALE; MOLARITY (1), MOLALITY (2),
C
C  THE FOLLOWING TABLE CONTAINS THE ALLOWABLE COMBINATIONS OF THESE FOUR
C   QUANTITIES.  A "-" INDICATES THAT THE VALUE IS IRRELEVANT.  FOR NASP,
C   THE VALUE "1" DENOTES ONE OR MORE SPECIES.
C
C------------------------------------------------------------------------
C             CONDITIONS            CALCN(S).REQD.       FLAGS
C
C       IAO   NASP  IEQ   ICN            NCAL        IAOC IOC  IAC
C------------------------------------------------------------------------     
C        0     -     0     -              2           0    1    0
C        1     -     -     -              1           1    0    0   
C        1     1     0     1              3, THEN 1   1    0    1
C        1     1     0     2              1           1    0    0
C        1     1     0     3              1           1    0    0
C        2     -     0     -              3           0    0    1
C------------------------------------------------------------------------
C
C  (IN THE CASE CALLING FOR TWO CALCULATIONS, THE FIRST IS NEEDED FOR
C   MATERIAL BALANCE CALCULATIONS WHEN THE MOLARITY SCALE IS USED.)
C
C
C                        MATERIAL BALANCE CORRECTION
C
C  A CORRECTION IS REQUIRED FOR THE CHANGE IN THE VOLUME OR THE WEIGHT OF
C   EACH PHASE THAT RESULTS WHEN MATERIAL IS TRANSFERRED FROM ONE PHASE TO
C   THE OTHER IN ORDER TO SATISFY MATERIAL BALANCE REQUIREMENTS. 
C
C  THE CORRECTION FACTORS USED FOR THIS PURPOSE, RMBO FOR THE NONAQUEOUS PHASE
C   AND RMBA FOR THE AQUEOUS PHASE, ARE EACH A RATIO OF THE AMOUNT (V0LUME,
C   WEIGHT, OR WEIGHT OF CONTAINED SOLVENT) OF A PHASE AFTER EXTRACTION TO
C   THAT BEFORE EXTRACTION.
C
C
C
      SUBROUTINE CALC(P,W,V,V0,SV,CS,G,
     &GO,GX0,VRT,PH,PHV,XF,GX,AC,X,XS,XCS,XCOS,XCAS,XC,	 
     &ROC,RAC,EOC,CEOC,REOC,SEOC,EAC,FOC,FAC,FOK,RFOK,SFOK,FAK,XCL,
     &FCL,EXC,WA,Y,B0,B1,B2,CP,ALP,BP,BG,BGP,CG,PCC,
     &ETHC,ETHCP,PAA,ETHA,ETHAP,PCCA,PCAA,DCA,PNC,PNA,PNN,AM,BM,CM,
     &ZA,ZFA,ZC,ZFC,CAP,COP,EP,PT,FX,PB,RZIR,SCPSI,KX,NZA,
     &NZPO,NCOP,NCAP,NIDT,NF,ICA,IAN,INU,IPT,IC,IDC,ICONV,IERROR,NPSI)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      EXTERNAL ENTFCN
      DIMENSION XCL(*),FCL(*),P(*),W(*),V(*),V0(*),SV(*),CS(*),
     &FOK(*),RFOK(5,*),SFOK(*),FAK(*),AM(*),BM(*),CM(*),B0(*),
     &B1(*),B2(*),CP(*),ALP(2,*),BP(*),BG(*),BGP(*),CG(*),
     &PCC(*),ETHC(*),ETHCP(*),PAA(*),ETHA(*),ETHAP(*),PCCA(*),
     &PCAA(*),DCA(*),PNC(*),PNA(*),PNN(*),ZA(*),ZFA(*),ZC(*),
     &ZFC(*),G(*),GO(*),GX0(*),VRT(*),PH(*),PHV(*),XF(*),
     &GX(*),AC(*),X(10,*),XC(*),CAP(*),COP(*),XS(11,*),XCS(*),
     &XCOS(*),XCAS(*),EXC(*),WA(*),ROC(8,*),RAC(8,*),EOC(8,*),
     &CEOC(*),REOC(5,8+NSO,*),SEOC(8+NSO,*),EAC(8,*),FOC(8,*),
     &FAC(8,*),EP(*),PT(8+NSO+NSA,*),FX(*),PB(*),RZIR(*),SCPSI(*)
      DIMENSION KX(*),NF(*),ICA(*),IAN(*),INU(*),NZA(*),NZPO(*),
     &NIDT(*),NCOP(8,*),NCAP(7,*),NPSI(*)
      COMMON/LGP/LP
      COMMON/MD/NSO,NCBO,IREF,NPH,NPA,NSA,NIO,NIA,NCA,NAN,NNU,NUP,NUM
      COMMON/DT/TC,NO,NOF,ICNR      
      COMMON/GRP/AFT,BFT,QFT,NCR(4),WS,VS,DS,GOS,WR,VR,GOR,RT
      COMMON/FIT/SIGMA,NC,NVP,LAM,LEV
      COMMON/WTR/P0,DW,DC,AP,RLKW      
      COMMON/EQ/DEN,AWG,AW,TWSA,TVSA,VWPA,VWPO,RMBA,RMBO
      COMMON/PRP/XFS,GXS,AS,RMCO,THSV,THS,TSAI,TSOI,TWSO,TVSO
      COMMON/RSDLS/SMINS
      COMMON/CYCLES/MXHCY,NTRY
      COMMON/ELP/EFT1,EFT2,ATM,SQI
      DATA FTOL,XTOL,TEN,CNT,THOU,DEXC,SMTL,TFCSM,NCSM/
     &1.0D-11,1.0D-11,10.0D+0,100.0D+0,1.0D+3,0.1D+0,1.0D-9,
     &1.0D+10,1000/
      DATA IT1,IT2,IT3,IT4/10,100,1000,10000/
      DATA ZERO,ONE,TWO/0.0D+0,1.0D+0,2.0D+0/
C
C  INITIALIZE FACTOR FOR ELECTROSTATIC CONTRIBUTION TO ACTIVITY
C   COEFFICIENTS AND SOLVENT ACTIVITY IN THE NONAQUEOUS PHASE
      EFT1=ONE
      EFT2=ONE
C
C  INITIALIZE CS ARRAY
      DO 3 I=1,(8+NSO+NSA)
      CS(I)=ZERO
3     CONTINUE
C
      IF(IDC.GT.0) GOTO 10
      IRES=NIDT(IPT)
C  DETERMINE TYPE OF DATA POINT
      IY=IRES/IT4
      IRES=IRES-IT4*IY
C  DETERMINE IF AQUEOUS PHASE IS PRESENT
      IAO=IRES/IT3
      IRES=IRES-IT3*IAO
C  DETERMINE IF INITIAL OR EQUILIBRIUM AQ.COMPOSITION IS SPECIFIED
      IEQ=IRES/IT2
      IRES=IRES-IT2*IEQ
C  DETERMINE IF INITIAL TOTAL ACIDITY, EQUIL.TOTAL ACIDITY, OR
C   ACTUAL H+ CONC.IS SPECIFIED
      IEH=IRES/IT1
C  DETERMINE CONCENTRATION SCALE
      ICN=IRES-IT1*IEH
C  DETERMINE THE NUMBER OF AQ.PROD.SPEC.NOT INVOLVING HA AND/OR B (NASP)
      NASP=0
      IF(NSA.GT.0) THEN
        DO 5 I=1,NSA
        IF(NCAP(3,I).EQ.0.AND.NCAP(4,I).EQ.0) THEN
C        THIS SPECIES DOES NOT CONTAIN NONAQ.REACTANTS A OR B
	  NASP=NASP+1
	END IF
5       CONTINUE
      END IF
C
C  DETERMINE THE EQUILIBRIUM CALCULATION(S) REQUIRED
      IAOC=0
      IOC=0
      IAC=0
C   ARE BOTH PHASES PRESENT?
      IF(IAO.EQ.1) IAOC=1
C   IS ONLY THE NONAQ.PHASE PRESENT?
      IF(IAO.EQ.0) IOC=1
C   IS ONLY THE AQ.PHASE PRESENT?     
      IF(IAO.EQ.2) IAC=1
C   ARE BOTH PHASES PRESENT, WITH AQ.PROD.SPEC.NOT INVOLVING HA OR B,
C    WITH INITIAL CONC.OF AQ.REACTANTS GIVEN, AND THE MOLARITY CONC.
C    SCALE SPECIFIED?     
      IF(IAO.EQ.1.AND.NASP.GT.0.AND.IEQ.EQ.0.AND.ICN.EQ.1) IAC=1
C  IF NEEDED, CALL SUBROUTINE TO PREPARE FOR NONAQ.ACT.COEF.CALCNS.
10    IF(IAOC.EQ.1.OR.IOC.EQ.1) CALL ORGAC(P,W,V,XF,GX,CS,G,
     &AC,GO,GX0,VRT,PH,PHV,NZPO,ICN)
C  IF NEEDED, CALL SUBROUTINE TO PREPARE FOR AQU.ACT.COEF.CALCULATIONS
      IF(IAOC.EQ.1.OR.IAC.EQ.1) CALL AQUAC(P,W,V,V0,SV,CS,AM,
     &BM,CM,B0,B1,B2,CP,ALP,BP,BG,BGP,CG,PCC,ETHC,ETHCP,
     &PAA,ETHA,ETHAP,PCCA,PCAA,DCA,PNC,PNA,PNN,ZA,ZFA,ZC,ZFC,G,
     &ICA,IAN,INU,NZA,ICN,NCAL)
C
C  BEGIN EQUILIBRIUM CALCULATION
      NCAL=0
      IF(IAOC.EQ.1.AND.IAC.EQ.0) THEN
C      A CALCULATION INVOLVING BOTH PHASES IS NEEDED
        NCAL=1
	GOTO 15
      END IF
      IF(IOC.EQ.1) THEN
C      A CALCULATION INVOLVING ONLY THE NONAQUEOUS PHASE IS NEEDED
        NCAL=2
	GOTO 15
      END IF
      IF(IAC.EQ.1) THEN
C      A CALCULATION INVOLVING THE INITIAL AQ.PHASE IS NEEDED
        NCAL=3
	GOTO 15
      END IF
      IF(NCAL.EQ.0) THEN
        WRITE(9,1000) IPT,NCAL
	WRITE(12,1000) IPT,NCAL
1000  FORMAT(2X,'PT.',I3,': NCAL  =',I2,': REQD.CALCNS. NOT ',
     &'DEFINED IN CALC'/)									
        PAUSE
	STOP
      END IF
C
15    CONTINUE
C
      IF(IAOC.EQ.1.AND.IDC.EQ.0) THEN
C      THIS IS THE FIRST ENTRY OF THIS CYCLE FOR A POINT WITH BOTH
C       PHASES PRESENT: ASSIGN INITIAL TOTAL SOLUTE CONTENTS
C       ASSIGN INITIAL ION CONCENTRATIONS IN THE AQUEOUS PHASE
        DO 20 I=3,6
        CS(I)=X(I,IPT)
20      CONTINUE
        CHI=NZA(4)*CS(5)+NZA(5)*CS(6)-NZA(1)*CS(3)+NZA(2)*CS(4)
	IF(CHI.GE.ZERO) THEN
	  CS(7)=CHI
	  CS(8)=ZERO
	ELSE
	  CS(7)=ZERO
	  CS(8)=-CHI
	END IF
        IF(NSA.GT.0) THEN
	  DO 25 I=1,NSA
	  CS(8+NSO+I)=ZERO
25        CONTINUE
        END IF
C       CALCULATE SOLUTE CONTENTS
        IF(ICN.EQ.1) THEN
C        CALCULATE INITIAL SOLUTE VOLUMES IN NONAQUEOUS PHASE
	  TSOI=V(1)*X(1,IPT)/NCR(1)+V(2)*X(2,IPT)/NCR(2)
C        SHOULD INITIAL APPARENT MOLAR VOLUMES BE CALCULATED HERE? 
          IF(IAC.EQ.0) THEN
C          CALCULATION WILL NOT BE MADE LATER. CALCULATE THEM HERE
	    CALL ENTAAC(P,W,V,V0,SV,CS,AM,
     &      BM,CM,B0,B1,B2,CP,ALP,BP,BG,BGP,CG,PCC,ETHC,ETHCP,
     &      PAA,ETHA,ETHAP,PCCA,PCAA,DCA,PNC,
     &      PNA,PNN,ZA,ZFA,ZC,ZFC,G,ICA,IAN,INU,NZA,ICN,NCAL)
            TSAI=TVSA
	  END IF
        ELSE
C        CALCULATE INITIAL SOLUTE WEIGHTS IN NONAQUEOUS PHASE
          TSOI=W(1)*X(1,IPT)/NCR(1)+W(2)*X(2,IPT)/NCR(2)
C        CALCULATE INITIAL SOLUTE WEIGHTS IN AQUEOUS PHASE
          TSAI=ZERO
          DO 30 I=3,8
          TSAI=TSAI+W(I)*CS(I)
30        CONTINUE
        END IF
      END IF
C
35    CALL EST(X,XS,XCS,XCOS,XCAS,XC,IPT, IDC,ICONV,NCAL,IST,IEQ,IEH)
C
C  CALL FCN TO DETERMINE NO. OF EQNS.TO BE SOLVED, AND TO MAKE
C   INITIAL ESTIMATE OF PRODUCT SPECIE CONCNS.
      CALL FCN(NEQ,XCL,FCL,IER,P,W,V,
     &V0,SV,CS,ROC,RAC,EOC,CEOC,REOC,SEOC,EAC,FOC,FAC,FOK,RFOK,
     &SFOK,FAK,AM,BM,CM,B0,B1,B2,CP,ALP,BP,BG,BGP,CG,
     &PCC,ETHC,ETHCP,PAA,ETHA,ETHAP,PCCA,PCAA,DCA,PNC,PNA,PNN,ZA,ZFA,
     &ZC,ZFC,G,GO,GX0,VRT,PH,PHV,XF,GX,AC,X,XC,CAP,
     &COP,RZIR,SCPSI,NCOP,NCAP,NF,ICA,IAN,INU,IPT,NZA,NZPO,IY,NCAL,
     &IEQ,IEH,ICN,IST,IDC,NPSI)
C
C     WRITE(17,4321) IREF,NEQ,(NF(K),K=1,NEQ)
C4321  FORMAT(/1X,'IREF,NEQ, NF:',2I5,10I3)
C
      IF(NEQ.EQ.0) THEN
        IF(NCAL.NE.3) 
     &   CALL ENTOAC(P,W,V,XF,GX,CS,G,AC,GO,GX0,VRT,PH,PHV,NZPO,ICN)
        IF(NCAL.NE.2)
     &   CALL ENTAAC(P,W,V,V0,SV,CS,AM,BM,CM,B0,B1,B2,CP,ALP,BP,
     &  BG,BGP,CG,PCC,ETHC,ETHCP,PAA,ETHA,ETHAP,PCCA,PCAA,DCA,PNC,
     &  PNA,PNN,ZA,ZFA,ZC,ZFC,G,ICA,IAN,INU,NZA,ICN,NCAL)     
        GOTO 95
      END IF
C
C   BEGIN REFINEMENT
C
      MAXCYC=0
      SMIN=SMINS/DBLE(NEQ)
      SMNSAV=SMIN
C    ASSIGN VALUES TO BE REFINED
      DO 40 I=1,NEQ
      XCL(I)=XC(NF(I))
40    CONTINUE
C
C      WRITE(*,1200)(XCL(I),I=1,NEQ)
C1200  FORMAT(1X,'    FROM EST: XCL = '/(1X,5F15.10))
C
C    DETERMINE IF SIMPLEX REFINEMENT OF XCL VALUES IS NEEDED
      CALL ENTFCN(NEQ,XCL,FCL,IER,P,W,V,
     &V0,SV,CS,ROC,RAC,EOC,CEOC,REOC,SEOC,EAC,FOC,FAC,FOK,RFOK,
     &SFOK,FAK,AM,BM,CM,B0,B1,B2,CP,ALP,BP,BG,BGP,CG,
     &PCC,ETHC,ETHCP,PAA,ETHA,ETHAP,PCCA,PCAA,DCA,PNC,PNA,PNN,ZA,ZFA,
     &ZC,ZFC,G,GO,GX0,VRT,PH,PHV,XF,GX,AC,X,XC,CAP,
     &COP,RZIR,SCPSI,NCOP,NCAP,NF,ICA,IAN,INU,IPT,NZA,NZPO,IY,NCAL,
     &IEQ,IEH,ICN,IST,IDC,NPSI)
      TFCL=ZERO
      DO 45 I=1,NEQ
      TFCL=TFCL+FCL(I)**2
      EXC(I)=DEXC
45    CONTINUE
      TFC=DSQRT(TFCL)/NEQ
C      WRITE(9,1049) NEQ,(NF(I),I=1,NEQ)
C      WRITE(17,1049) NEQ,(NF(I),I=1,NEQ)
C1049  FORMAT(1X,'NEQ =',I3/1X,' NF = ',8I9)      
C      WRITE(9,1070) (XCL(I),I=1,NEQ)
C      WRITE(17,1070) (XCL(I),I=1,NEQ)
C1070  FORMAT(1X,'XCL = ',1P8E9.1)
C      WRITE(9,1050) (FCL(I),I=1,NEQ)
C      WRITE(17,1050) (FCL(I),I=1,NEQ)
C1050  FORMAT(1X,'FCL = ',1P8E9.1)      
48    WRITE(17,1001) IPT,NCAL,SMIN,TFC
C      WRITE(9,1001) IPT,NCAL,SMIN,TFC
1001  FORMAT(2X,'PT.',I3,': NCAL  =',I2,
     &':  DESIRED RES. =',F7.4,': INITIAL RES. =',F9.3)
C
C      IF((TFC/SMIN).GT.TEN) THEN
C        WRITE(17,2113)(NF(I),I=1,NEQ)
C2113    FORMAT(1X,'NF = ',(10I4))
C        WRITE(17,2120)(XCL(I),I=1,NEQ)
C2120    FORMAT(1X,'XCL = '/(1X,5F15.10))
C        WRITE(17,2121) (FCL(I),I=1,NEQ)
C2121    FORMAT(1X,'FCL = ',1P8E9.2)
C      END IF
C     
      IF(TFC.LE.SMIN) GO TO 80
C  REFINE VALUES BY USE OF THE SIMPLEX PROCEDURE
      NEQU=NEQ
      SMTOL=SMTL
      NCYSM=NCSM
      CALL SMPLX(NEQU,EXC,SMTOL,NCYSM,EP,PT,FX,KX,PB)
      NSX=0
      TFCS=TFCSM
50    DO 70 KK=1,IT2
      IF(NSX.EQ.0) GOTO 60
      CALL ENTFCN(NEQ,XCL,FCL,IER,P,W,V,
     &V0,SV,CS,ROC,RAC,EOC,CEOC,REOC,SEOC,EAC,FOC,FAC,FOK,RFOK,
     &SFOK,FAK,AM,BM,CM,B0,B1,B2,CP,ALP,BP,BG,BGP,CG,
     &PCC,ETHC,ETHCP,PAA,ETHA,ETHAP,PCCA,PCAA,DCA,PNC,PNA,PNN,ZA,ZFA,
     &ZC,ZFC,G,GO,GX0,VRT,PH,PHV,XF,GX,AC,X,XC,CAP,
     &COP,RZIR,SCPSI,NCOP,NCAP,NF,ICA,IAN,INU,IPT,NZA,NZPO,IY,NCAL,
     &IEQ,IEH,ICN,IST,IDC,NPSI)
      TFCL=ZERO
      DO 55 I=1,NEQ
      TFCL=TFCL+FCL(I)**2
55    CONTINUE
      TFC=DSQRT(TFCL)/NEQ
60    NSX=NSX+1
      IF(TFC.GE.TFCS) GOTO 65
      TFCS=TFC
65    CALL SMPLX(0,XCL,TFC,IRX,EP,PT,FX,KX,PB)
      IF(IRX.NE.0.OR.TFC.LE.SMIN) GOTO 75
70    CONTINUE
      WRITE(9,1002) IPT,NSX,TFC
      GOTO 50
75    WRITE(17,1002) IPT,NSX,TFC 
      WRITE(9,1002) IPT,NSX,TFC 
1002  FORMAT(3X,'POINT NO.',I3,': SMPLX CYCLES  = ',I4,
     &': RESIDUAL. =',F8.3)
C   PERFORM FINAL REFINEMENT OF CONCENTRATIONS
80    MAXFEV=MXHCY
      LWA=NEQ*(3*NEQ+7)/2
      IERROR=0
C
      CALL HYBRD2(NEQ,ENTFCN,XCL,FCL,FTOL,XTOL,MAXFEV,
     &IER,LWA,WA,P,W,V,V0,SV,CS,ROC,RAC,EOC,CEOC,REOC,
     &SEOC,EAC,FOC,FAC,FOK,RFOK,SFOK,FAK,AM,BM,CM,B0,B1,B2,
     &CP,ALP,BP,BG,BGP,CG,PCC,ETHC,ETHCP,PAA,ETHA,ETHAP,PCCA,PCAA,
     &DCA,PNC,PNA,PNN,ZA,ZFA,ZC,ZFC,G,GO,GX0,VRT,PH,PHV,
     &XF,GX,AC,X,XC,CAP,COP,RZIR,SCPSI,NCOP,NCAP,NF,ICA,IAN,
     &INU,IPT,NZA,NZPO,IY,NCAL,IEQ,IEH,ICN,IST,IDC,NPSI)
      MAXCYC=MAXCYC+MAXFEV
      IF(IER.LE.3) GO TO 85
      IF(IER.EQ.4.OR.IER.EQ.5.OR.IER.EQ.6) THEN
        TFCL=ZERO
        DO 83 I=1,NEQ
        TFCL=TFCL+FCL(I)**2
        EXC(I)=DEXC
83      CONTINUE
        TFC=DSQRT(TFCL)/NEQ
        IF(SMIN.GT.(SMNSAV/TEN)) SMIN=SMIN/TWO
        GOTO 48
      END IF
      IF(NCAL.EQ.1) WRITE(9,1003) IPT,IER
      IF(NCAL.EQ.1) WRITE(12,1003) IPT,IER
      IF(NCAL.EQ.2) WRITE(9,1004) IPT,IER
      IF(NCAL.EQ.2) WRITE(12,1004) IPT,IER
      IF(NCAL.EQ.3) WRITE(9,1005) IPT,IER
      IF(NCAL.EQ.3) WRITE(12,1005) IPT,IER
1003  FORMAT(1X,' BOTH PHASES PRESENT: COMPUTATION DID NOT CONVERGE',
     &' FOR POINT',I3/5X,'ERROR CODE = ',I2/)
1004  FORMAT(1X,' ORG.PHASE ONLY: COMPUTATION DID NOT CONVERGE',
     &' FOR POINT',I3/5X,'ERROR CODE = ',I2/)
1005  FORMAT(1X,' AQ.PHASE ONLY: COMPUTATION DID NOT CONVERGE',
     &' FOR POINT',I3/5X,'ERROR CODE = ',I2/)
C      WRITE(9,1006) IPT, NIDT(IPT),(X(KK,IPT),KK=1,10)
C1006  FORMAT(/1X,'FOR POINT',I3,' NIDT =',I7/
C     &1X,'X(1) - X(10):'/1X,1P5E15.5/1X,5E15.5)
C      WRITE(9,1007) (XCL(KK),KK=1,NEQ)
C1007  FORMAT(1X,'XCL(1) - XCL(NEQ):'/
C     &1X,1P5E15.5/1X,5E15.5/1X,5E15.5)
C      WRITE(9,1008) (FCL(KK),KK=1,NEQ)
C1008  FORMAT(1X,'FCL(1) - FCL(NEQ):'/
C     &1X,1P5E15.5/1X,5E15.5/1X,5E15.5)
      IERROR=1
      PAUSE
      STOP
C  CALCULATE EQUILIBRIUM SPECIE CONCENTRATION
85    CONTINUE
C
C      WRITE(*,1210)(NF(I),I=1,NEQ)
C1210  FORMAT(1X,'NF = ',(10I4))
C      WRITE(*,1220)(XCL(I),I=1,NEQ)
C1220  FORMAT(1X,'AFTER HYBRD2: XCL = '/(1X,5F15.10))
C      PAUSE
C           
      CALL ENTFCN(NEQ,XCL,FCL,IER,P,W,V,
     &V0,SV,CS,ROC,RAC,EOC,CEOC,REOC,SEOC,EAC,FOC,FAC,FOK,RFOK,
     &SFOK,FAK,AM,BM,CM,B0,B1,B2,CP,ALP,BP,BG,BGP,CG,
     &PCC,ETHC,ETHCP,PAA,ETHA,ETHAP,PCCA,PCAA,DCA,PNC,PNA,PNN,ZA,ZFA,
     &ZC,ZFC,G,GO,GX0,VRT,PH,PHV,XF,GX,AC,X,XC,CAP,
     &COP,RZIR,SCPSI,NCOP,NCAP,NF,ICA,IAN,INU,IPT,NZA,NZPO,IY,NCAL,
     &IEQ,IEH,ICN,IST,IDC,NPSI)
      DO 90 I =1,NEQ
      XC(NF(I))=XCL(I)
90    CONTINUE
95    CONTINUE
C  SAVE SPECIE CONCNS.
      IF(NCAL.EQ.1) THEN
        K=(8+NSO+NSA)*(IPT-1)
        DO 100 I=1,(8+NSO+NSA)
	K=K+1
        XCS(K)=XC(I)
100     CONTINUE
      END IF
      IF(NCAL.EQ.2) THEN
        K=(8+NSO+NSA)*(IPT-1)
        DO 105 I=1,(8+NSO+NSA)
	K=K+1
        XCOS(K)=XC(I)
105     CONTINUE
      END IF
      IF(NCAL.EQ.3) THEN
        K=(6+NSA)*(IPT-1)
        DO 110 I=3,8
	K=K+1
        XCAS(K)=XC(I)
110     CONTINUE
        IF(NSA.GT.0) THEN
          DO 115 I=(9+NSO),(8+NSO+NSA)
          K=K+1
          XCAS(K)=XC(I)
115       CONTINUE
        END IF
      END IF
C      WRITE(9,1010) MAXFEV
      WRITE(17,1010) MAXFEV
1010  FORMAT(3X,'HYBRD2 CYCLES =',I4)
      IF(NCAL.EQ.3.AND.IAOC.EQ.1) THEN
C      CALCULATION COMPLETE FOR INITIAL AQ.PHASE,
C       PERFORM EQ.CALCNS.WITH BOTH PHASES PRESENT
        NCAL=1
	GOTO 35
      END IF
      IF(IC.LT.(NC+1)) GOTO 120
C
C      WRITE(16,1015)
C1015  FORMAT(/)
C
      IF(IAO.EQ.1) THEN
C      BOTH PHASES ARE PRESENT
        IF(ICN.EQ.1) THEN
          WRITE(16,1020) IPT,ICN,DEN,TSAI,TVSA,VWPA,RMBA,
     &    TSOI,TVSO,VWPO,RMBO
1020      FORMAT(1X,I3,2X,I2,2X,F11.6,3(2X,F9.3),2X,F11.6/
     &    21X,3(2X,F9.3),2X,F11.6)
        ELSE
          WRITE(16,1021) IPT,ICN,TSAI,TWSA,VWPA,RMBA,
     &    TSOI,TWSO,VWPO,RMBO
1021      FORMAT(1X,I3,2X,I2,13X,3(2X,F9.3),2X,F11.6/
     &    21X,3(2X,F9.3),2X,F11.6)
        END IF
      END IF
      WRITE(14,1023) IPT,(CS(I),I=1,(8+NSO+NSA))
1023  FORMAT(1X,I3,1P5E14.7/4X,3E14.7/4X,5E14.7/4X,5E14.7/
     &4X,5E14.7/4X,5E14.7/4X,5E14.7/4X,5E14.7)
      IF((NSO+NSA).EQ.0.OR.(NSO+NSA).EQ.5.OR.(NSO+NSA).EQ.10.OR.
     &(NSO+NSA).EQ.15.OR.(NSO+NSA).EQ.20) BACKSPACE 14
      WRITE(15,1024) IPT,(G(I),I=1,(8+NSO+NSA)),EFT1
1024  FORMAT(1X,I3,0P8F9.5/4X,8F9.5/4X,8F9.5/4X,8F9.5/4X,8F9.5)
      IF((NSO+NSA).EQ.0.OR.(NSO+NSA).EQ.8.OR.(NSO+NSA).EQ.16)
     &BACKSPACE 15
      WRITE(15,1025) (AC(I),I=1,(2+NSO))
1025  FORMAT(4X,1P8E9.2/4X,8E9.2/4X,8E9.2/4X,8E9.2/4X,8E9.2)
      IF(NSO.EQ.6.OR.NSO.EQ.14.OR.NSO.EQ.22) BACKSPACE 15
      WRITE(15,1026) AW,AS,EFT2
1026  FORMAT(4X,3F9.6)  
120   CONTINUE
C      WRITE(15,1100) (XF(I),I=1,(2+NSO))
C1100  FORMAT(4X,1P5E13.5/(4X,1P5E13.5))
C      WRITE(15,1101) (GX(I),I=1,(2+NSO))
C1101  FORMAT(4X,0P8F9.3/(4X,8F9.3))
C
C  CALCUALTE AN EXTRACTION COEFFICIENT OF M
      IF(IY.EQ.1.OR.IY.EQ.3.OR.IY.EQ.5.OR.IY.EQ.7)
     &CALL DCALC(IY,CS,NCAP,NCOP,Y)
C  CALCUALTE A TOTAL CONC. IN NONAQUEOUS PHASE
      IF(IY.EQ.2.OR.IY.EQ.4.OR.IY.EQ.6.OR.IY.EQ.8)
     &CALL TCALC(IY,CS,NCOP,Y)      
C  CALCULATE THE WATER CONTENT OF THE NONAQUEOUS PHASE
      IF(IY.EQ.10) CALL WCALC(CS,NCOP,Y)
C  CALCUALTE LOG EXTRACTION COEFFICIENT OF M
      IF(IY.EQ.11.OR.IY.EQ.13.OR.IY.EQ.15.OR.IY.EQ.17)
     &CALL LDCAL(IY,CS,NCAP,NCOP,Y)
C  CALCUALTE LOG TOTAL CONC. IN NONAQUEOUS PHASE
      IF(IY.EQ.12.OR.IY.EQ.14.OR.IY.EQ.16.OR.IY.EQ.18)
     &CALL LTCAL(IY,CS,NCOP,Y)      
C  CALCULATE THE PARTICLE CONCENTRATION
      IF(IY.EQ.20) CALL PCALC(ICN,Y)
C  CALCULATE THE HEAT OF MIXING OF TWO NONAQUEOUS SOLUTIONS
      IF(IY.EQ.30) CALL HCALC(IPT,P,X,CS,G,XF,ICN,Y)
C  CALCULATE THE SPECTRAL ABSORBENCE 
      IF(IY.EQ.40.OR.IY.EQ.41) CALL ECALC(IPT,P,X,CS,Y)
C  CALCULATE THE ACTIVITY COEFFICIENT OF A COMPONENT, THE OSMOTIC 
C   COEFFICIENT, OR THE WATER ACTIVITY OF THE AQUEOUS PHASE
      IF(IY.GE.60.AND.IY.LT.70) CALL GCALC(IY,IPT,NZA,X,CS,G,Y)
C  CALCULATE A QUANTITY FROM A USER SUBROUTINE
      IF(IY.GE.70.AND.IY.LT.80)
     &CALL USRCAL(IY,IPT,ICN,P,CS,X,G,XF,NCAP,NCOP,NZA,Y,V0,SV)
      RETURN
      END
C
C
C
C  SUBROUTINE TO CALCULATE NONAQ. ACTIVITY COEFFICIENTS ACCORDING
C   TO THE TREATMENT OF HILDEBRAND AND SCOTT, INCLUDING THE FLORY-
C   HUGGINS CORRECTION RECOMMENDED BY HENLEY AND SEADER (1981)
C
C   VOLUME ADDITIVITY OF ALL SPECIES IS ASSUMED
C
C
C                            GLOSSARY
C
C      AC   =  ACTIVITY OF NONAQ. SOLUTE (X SCALE)
C      AS   =  ACTIVITY OF NONAQ. SOLVENT
C      ATM  =  DEBEYE-HUCKEL TERM (IN NUMERATOR)
C      BTM  =  DEBEYE-HUCKEL TERM (IN DENOMINATOR)
C      CI   =  IONIC STRENGTH (MOLARITY SCALE)
C      CS   =  CURRENT CONCENTRATION OF A SOLUTE
C      DS   =  DIELECTRIC CONSTANT OF THE NONAQUEOUS SOLVENT
C      EFT  =  ELECTROSTATIC TERM IN ACTIVITY COEFFICIENT OF AN ION
C      EFT1 =  ELECTROSTATIC TERM IN ACTIVITY COEFFICIENT OF A
C      EFT2 =  ELECTROSTATIC TERM IN ACTIVITY SOLVENT
C      FHR  =  FLORY-HUGGINS RATIO (V/(SUM(XV))
C      G    =  ACTIVITY COEFFICIENT OF A SOLUTE ON THE GIVEN
C               CONCENTRATION SCALE
C      GAV  =  AVERAGE SOLUBILITY PARAMETER FOR MIXT., WEIGHTED BY
C               VOLUME FRACTION
C      GO   =  SOLUBILITY PARAMETER OF SOLUTE (CAL**1/2/CM**3/2)  
C      GOR  =  SOLUBILITY PARAMETER OF THE REFERENCE SOLUTE
C      GOS  =  SOLUBILITY PARAMETER OF SOLVENT
C      GTM  =  DISTANCE TERM IN BJERRUM MODEL
C      GX   =  MOLE FRACTION ACT.COEF.OF SOLUTE
C      GXS  =  MOLE FRACTION ACT.COEF.OF SOLVENT
C      GX0  =  MOLE FRACTION ACT.COEF.OF SOLUTE AT INFINITE DILUTION
C      ICN  =  CONCENTRATION SCALE INDICATOR: 1 = MOLARITY, 2 = 
C               MOLALITY, 3 = MOL / KG SOLN.
C      IREF =  SPECIE NO.ASSIGNED TO IONIC NONAQUEOUS PRODUCT SPECIES
C               THAT SERVES AS THE REFERENCE ION IN FORMATION REACTIONS.
C               IF IREF IS ZERO, THEN NO NONAQUEOUS IONS ARE PRESENT.
C      NIO  =  INDICATOR OF CALCULATED (0) OR ASSUMED UNIT (1) VALUES
C               FOR NONAQ.ACT.COEFFICIENTS, G
C      PH   =  VOLUME FRACTION OF A SOLUTE 
C      PHS  =  VOLUME FRACTION OF SOLVENT
C      RMCO =  RATIO OF SOLUTE CONCN. IN MOLALITY TO THE GIVEN CONCN.
C      RXC  =  RATIO OF MOLE FRACTION TO GIVEN CONCENTRATION
C      RXCF =  LIMIT OF MOLE FRACTION TO GIVEN CONC. AT INFINITE
C               DILUTION
C      SPV  =  TOTAL MOLES OF SPECIES (INCLUDING SOLVENT) PER CM3 SOLN.
C      T    =  TEMPERATURE (K)
C      TC   =  TEMPERATURE (C)
C      TCS  =  TOTAL CONCENTRATION OF SOLUTES (MOLES PER LITER SOLN,
C               PER KG OF SOLVENT, OR PER KG OF SOLUTION)
C      THS  =  HEAT TERM DUE TO NON-IDEAL EFFECTS FROM THE SOLUTES
C               (PER G SOLVENT, CALCULATED IN SUB.ORGAC)
C      THSV =  HEAT TERM DUE TO NON-IDEAL EFFECTS FROM THE SOLVENT
C               (PER G SOLVENT, CALCULATED IN SUB.ORGAC)
C      TVL  =  VOLUME OF SOLUTION (CM3 PER LITER OF SOLUTION, PER KG
C               OF SOLVENT, OR PER KG OF SOLUTION) 
C      TVSO =  TOTAL VOLUME OF SOLUTES (CM3 PER LITER SOLN,
C               PER KG OF SOLVENT, OR PER KG OF SOLUTION)
C      TWL  =  WEIGHT OF SOLUTION (G PER LITER OF SOLUTION, PER KG
C               OF SOLVENT, OR PER KG OF SOLUTION) 
C      TWSO =  TOTAL WEIGHT OF SOLUTES (G PER LITER SLN.,PER KG OF
C               SOLVENT, OR PER KG OF SOLUTION)
C      V    =  MOLAR VOLUME OF A SOLUTE
C      VR   =  MOLAR VOLUME OF THE REFERENCE SOLUTE (CM3/MOLE)
C      VS   =  MOLAR VOLUME OF THE NONAQ.SOLVENT (CM3/MOLE)
C      W    =  MOLECULAR WT. OF A SOLUTE      
C      WR   =  MOLECULAR WT. OF THE REFERENCE SOLUTE
C      WS   =  MOLECULAR WT. OF THE NONAQ.SOLVENT
C      XF   =  MOLE FRACTION OF A SOLUTE
C      XFS  =  MOLE FRACTION OF SOLVENT
C      ZPO  =  CHARGE ON A PRODUCT SPECIES (WITHOUT SIGN)
C
C
C
      SUBROUTINE ORGAC(P,W,V,XF,GX,CS,G,AC,GO,GX0,VRT,PH,PHV,
     &NZPO,ICN)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      COMMON/MD/NSO,NCBO,IREF,NPH,NPA,NSA,NIO,NIA,NCA,NAN,NNU,NUP,NUM
      COMMON/DT/TC,NO,NOF,ICNR      
      COMMON/GRP/AFT,BFT,QFT,NCR(4),WS,VS,DS,GOS,WR,VR,GOR,RT
      COMMON/PRP/XFS,GXS,AS,RMCO,THSV,THS,TSAI,TSOI,TWSO,TVSO
      COMMON/ELP/EFT1,EFT2,ATM,SQI
      DIMENSION P(*),W(*),V(*),XF(*),GX(*),CS(*),G(*),AC(*),GO(*),
     &GX0(*),VRT(*),PH(*),PHV(*),NZPO(*)
      DATA ZERO,ONE,TWO,THREE,EIGHT,THOU/
     &0.0D+00,1.0D+0,2.0D+0,3.0D+0,8.0D+0,1.0D+3/
      DATA TK,AVN,BK,E,PI/
     &273.15D+0,6.022045D+23,1.380662D-16,4.803242D-10,3.1415926D+0/
      VSRT=VS/RT
C  ASSIGN SOLUBILITY PARAMETER VALUES
      K=NSO*(1+NPH+NPA)+2*(NPH+NPA)+NSA
      DO 10 I=1,(NSO+2)
      K=K+1
      IF(NIO.NE.1) GO(I)=P(K)
      IF(NIO.EQ.1) GO(I)=GOS
10    CONTINUE
C  EVALUATE CONSTANTS IN EXPRESSION FOR ACTIVITY COEFFICIENTS OF
C   IONIC SPECIES, IF NEEDED
      IF(IREF.NE.0.AND.NIO.EQ.0) THEN
C      ACTIVITY COEFS.WILL BE CALCULATED FOR IONIC SPECIES
        T=TC+TK
	ATM=DSQRT(WS/VS)*AFT/(DS*T)**(THREE/TWO)
	BTM=DSQRT(WS/VS)*BFT/DSQRT(DS*T)
	QTM=QFT/(DS*T)
C      CALCULATED SMALLEST PRODUCT OF Z+ AND Z-
        NZPOC=10
        NZPOA=-10
	DO 12 I=1,NSO
	IF(NZPO(I).GT.0.AND.NZPO(I).LT.NZPOC) NZPOC=NZPO(I)
        IF(NZPO(I).LT.0.AND.NZPO(I).GT.NZPOA) NZPOA=NZPO(I)
12      CONTINUE
        ZPOCA=DABS(DBLE(NZPOC)*DBLE(NZPOA))
C
C        WRITE(9,4321) ATM,BTM,QTM,ZPOCA
C4321    FORMAT(/1X,'            ATM            BTM            QTM',
C     &  '   ZPOCA'/1X,1P3E15.8,0PF8.1)
C        PAUSE
C      CALCULATE CONSTANT TERM ASSOCIATED WITH CONTRIBUTION OF IONS TO
C       SOLVENT ACTIVITY
        SEFCT=TWO*WS*ATM/(THOU*(ZPOCA*ATM)**3)
      END IF
C  CALCULATE  MOLE FRACTION ACTIVITY COEFFICIENT AT INFINITE DILUTION
C   FOR EACH SOLUTE SPECIES AND OTHER QUANTITIES THAT MAY BE NEEDED
      DO 15 I=1,(2+NSO)
      IF(I.LE.2) K=I
      IF(I.GT.2) K=6+I
      VRT(I)=V(K)/RT
      FHR=V(K)/VS
      GX0(I)=FHR*DEXP(VRT(I)*(GO(I)-GOS)**2+ONE-FHR)
15    CONTINUE
      RETURN
C
C  ENTRY POINT FOR CALCULATING CURRENT ACTIVITY COEFFICIENTS
      ENTRY ENTOAC(P,W,V,XF,GX,CS,G,AC,GO,GX0,VRT,PH,PHV,NZPO,ICN)
C
C      WRITE(*,1003) CS(1),CS(2),(CS(KK),KK=9,(8+NSO))
C1003  FORMAT(1X,'CS:'/(1X,1P6E13.2/))
C
C      PAUSE
C
C  CALCUALTE THE VOLUME FRACTION OF EACH SOLUTE SPECIES, THE AVERAGE
C   SOLUBILITY PARAMETER, AND THE TOTAL SOLUTE CONCENTRATION
      TCS=ZERO
      TVSO=ZERO
      TWSO=ZERO
      CI=ZERO
      DO 20 I=1,(2+NSO)
      IF(I.LE.2) K=I
      IF(I.GT.2) K=6+I
      TCS=TCS+CS(K)
      IF(I.GT.2.AND.NZPO(I-2).NE.0) THEN
C      ADD CONTRIBUTION TO IONIC STRENGTH
        CI=CI+(DBLE(NZPO(I-2)))**2*CS(K)/TWO
      END IF
      TVSO=TVSO+V(K)*CS(K)
      TWSO=TWSO+W(K)*CS(K)
20    CONTINUE
C
C      WRITE(9,3210) TVSO,TWSO,(CS(K),K=9,(8+NSO))
C      WRITE(12,3210) TVSO,TWSO,(CS(K),K=9,(8+NSO))
C3210  FORMAT(/1X,'        TVSO        TWSO       CS(I)'/
C     &1X,1P6E12.3/(25X,4E12.3/))
C      PAUSE
C
      IF(ICN.EQ.1) THEN
        TVL=THOU
	TWL=(THOU-TVSO)*WS/VS+TWSO
	RMCO=TVL/(TWL-TWSO)
      END IF
      IF(ICN.EQ.2) THEN
        TVL=THOU*VS/WS+TVSO
	TWL=THOU+TWSO
	RMCO=ONE
      END IF
      IF(ICN.EQ.3) THEN
        TVL=(THOU-TWSO)*VS/WS+TVSO
	TWL=THOU
	RMCO=TWL/(TWL-TWSO)
      END IF
C
      IF(CI.GT.ZERO) THEN
C      IONS ARE PRESENT: CORRECT CI TO MOLALITY SCALE IF NECESSARY
C       AND CALCULATE THE DEBEYE-HUCKEL FACTOR
        CI=CI*RMCO
        SQI=DSQRT(CI)
	DNTM=ONE+ZPOCA*BTM*QTM*SQI
        DHFT=-ATM*SQI/DNTM
	EFT1=DEXP(DHFT)
      END IF
C
      PHS=ONE
      GAV=ZERO
      DO 25 I=1,(2+NSO)
      IF(I.LE.2) K=I
      IF(I.GT.2) K=6+I
      PH(I)=V(K)*CS(K)/TVL
C
C      WRITE(*,1005) I,K,TVL,V(K),CS(K),PH(I)
C      PAUSE
C
      IF(PH(I).GT.ONE) THEN
C      A VOLUME FRACTION CAN NOT EXCEED UNITY
        WRITE(*,1005) I,K,TVL,V(K),CS(K),PH(I)
1005    FORMAT(/1X'A VOLUME FRACTION HAS EXCEEDED UNITY'/
     &  1X,'  I  K          TVL      V(K)        CS(K)        PH(I)'/
     &  1X,2I3,F10.2,1P3E13.2/)
        PAUSE
        STOP
      END IF
      GAV=GAV+PH(I)*GO(I)
      PHS=PHS-PH(I)
      IF(PHS.LT.ZERO) THEN
C      THE VOLUME FRACTION OF THE SOLVENT CAN NOT BE NEGATIVE
        WRITE(9,1006) I,PHS
 1006   FORMAT(/1X,'THE CALCULATED VOLUME FRACTION OF THE SOLVENT IS',
     &   F10.3)
        PAUSE
        STOP
      END IF
25    CONTINUE
      GAV=GAV+PHS*GOS
C  CALCULATE MOLE FRACTIONS
      SPV=ZERO
      DO 30 I=1,(2+NSO)
      IF(I.LE.2) K=I
      IF(I.GT.2) K=6+I
      PHV(I)=PH(I)/V(K)
      SPV=SPV+PHV(I)
30    CONTINUE
      SPV=SPV+PHS/VS
      XFS=ONE
      DO 35 I=1,(2+NSO)
      IF(I.LE.2) K=I
      IF(I.GE.3) K=6+I
      XF(I)=PHV(I)/SPV
      XFS=XFS-XF(I)
35    CONTINUE
C  ASSIGN RATIO (X/C), THE LIMIT OF (X/C) AT INFINITE DILUTION,
C   AND THE RATIO OF THESE TWO QUANTITIES
      IF(ICN.EQ.1) THEN
        RXC=ONE/(TCS+(THOU-TVSO)/VS)
	RXC0=VS/THOU
      END IF
      IF(ICN.EQ.2) THEN
        RXC=ONE/(TCS+THOU/WS)
	RXC0=WS/THOU
      END IF
      IF(ICN.EQ.3) THEN
        RXC=ONE/(TCS+(THOU-TWSO)/WS)
	RXC0=WS/THOU
      END IF
      RXCF=RXC/RXC0
      IF(NIO.NE.1) THEN
C      A NON-IDEAL SOLUTION IS ASSUMED
	HS=ZERO
        DO 40 I=1,(2+NSO)
        IF(I.LE.2) K=I
        IF(I.GE.3) K=6+I
C      CALCULATE SOLUTE ACTIVITIES ON MOLE FRACTION SCALE
        FHR=V(K)*SPV
        GX(I)=FHR*DEXP(VRT(I)*(GO(I)-GAV)**2+ONE-FHR)
        IF(NZPO(I-2).EQ.0) AC(I)=XF(I)*GX(I)
        IF(NZPO(I-2).NE.0) AC(I)=ZERO
C      CALCULATE SOLUTE ACTIVITY COEFFICIENTS ON GIVEN CONC. SCALE
C        WRITE(*,1111) I,K,GX(I),GX0(I)
C1111  FORMAT(1X,'IN ORGAC: I,K,GX(I),GX0(I):',2I4,1P2E12.3/)
        G(K)=RXCF*GX(I)/GX0(I)
C
        IF(I.GT.2.AND.NZPO(I-2).NE.0.AND.SQI.GT.ZERO) THEN
C        THIS IS AN IONIC SPECIES: ADD ELECTROSTATIC FACTOR TO G(K)
          EFT=DEXP(DHFT*(DBLE(NZPO(I-2)))**2)
          G(K)=G(K)*EFT
C          WRITE(9,4322) I,NZPO(I-2),SQI,ZPOCA,EFT
C4322      FORMAT(/1X,'  I NZPO(I-2)      SQI    ZPOCA      EFT'/
C     &    1X,I3,I10,3F9.5/)
C          PAUSE
        END IF
C
C      CALCULATE HEAT TERMS ARISING FROM SOLUTE NON-IDEALITY, 
C       CONVERTING CONC. TO MOLALITY IF NECESSARY
	HS=HS+RMCO*CS(K)*V(K)*((GO(I)-GAV)**2-(GO(I)-GOS)**2)
40      CONTINUE
	THS=HS/THOU
C      CALCULATE SOLVENT ACTIVITY
        FHR=VS*SPV
	DGS=(GOS-GAV)**2	
        GXS=FHR*DEXP(VSRT*DGS+ONE-FHR)
        AS=XFS*GXS	
	IF(CI.GT.ZERO) THEN
C        CALCULATE EFFECT OF IONIC STRENGTH ON SOLVENT ACTIVITY
	  EFT2=DEXP(SEFCT*(DNTM-ONE/DNTM-TWO*DLOG(DNTM)))
C
C	   WRITE(*,2001) EFT2
C2001      FORMAT(1X,'ELEXT.CORR.FACTOR FOR SOLVENT =',F10.7)
C          PAUSE
C
          AS=AS*EFT2
        END IF
C      CALCULATE HEAT TERM ARISING FROM THE SOLVENT
        THSV=(VS/WS)*DGS
      END IF
      IF(NIO.EQ.1) THEN
C      AN IDEAL SOLUTION IS ASSUMED      
        DO 50 I=1,(2+NSO)
        IF(I.LE.2) K=I
        IF(I.GT.2) K=6+I
C      ASSIGN UNIT SOLUTE ACTIVITY COEFFICIENTS
        G(K)=ONE
C      CALCULATE SOLUTE ACTIVITIES ON MOLE FRACTION SCALE
        IF(ICN.EQ.1) AC(I)=PH(I)
        IF(ICN.EQ.2) AC(I)=ZERO
        IF(ICN.EQ.3) AC(I)=CS(K)*W(K)/THOU
50      CONTINUE
C      CALCULATE SOLVENT ACTIVITY
        IF(ICN.EQ.1) THEN
          AS=DEXP((VS*TCS/(THOU*(ONE-PHS)))*DLOG(PHS))
          RETURN
        END IF
        IF(ICN.EQ.2) THEN
          AS=DEXP(-(WS/THOU)*TCS)
          RETURN
        END IF
        IF(ICN.EQ.3) THEN
          AS=DEXP((WS*TCS/TWSO)*DLOG(ONE-TWSO/THOU))
          RETURN
        END IF
C      ASSIGN ZERO TO HEAT TERMS
        THS=ZERO
	THSV=ZERO
      END IF
C      WRITE(12,1113) CS(1),CS(2),(CS(K),K=9,(NSO+8))
C1113  FORMAT(1X,'    CS (HA B NSO) ='/1X,7F10.6)
C      WRITE(12,1114) PH(1),PH(2),(PH(K),K=3,(NSO+2))
C1114  FORMAT(1X,'    PH (HA B NSO) ='/1X,7F10.6)
C      WRITE(12,1115) XF(1),XF(2),(XF(K),K=3,(NSO+2))
C1115  FORMAT(1X,'    XF (HA B NSO) ='/1X,7F10.6)
C      WRITE(12,1116) GO(1),GO(2),(GO(K),K=3,(NSO+2))
C1116  FORMAT(1X,'    GO (HA B NSO) ='/1X,7F10.3)
C      WRITE(12,1117) G(1),G(2),(G(K),K=9,(NSO+8))
C1117  FORMAT(1X,'     G (HA B NSO) ='/1X,7F10.6)
C      WRITE(12,1118) VS,WS,PHS,XFS,GOS,GAV,VSRT,AS
C1118  FORMAT(1X,'VS,WS,PHS,XFS,GOS,GAV,VSRT,AS ='/
C    & 1X,2F9.3,2F9.5,2F9.3,2F9.5/)
C      PAUSE
      RETURN
      END
C
C
C
C  SUBROUTINE TO CALCULATE THE OSMOTIC COEFFICIENT AND ACTIVITY
C   COEFFICIENTS OF SPECIES IN AN ELECTROLYTE MIXTURE FROM THE
C   TREATMENT OF PITZER
C
C
C                             GLOSSARY
C
C     ALP         =  PITZER CONSTANTS ALPHA
C     AM          =  ANION CONCENTRATION
C     AP          =  DEBYE-HUCKEL CONSTANT A(PHI)
C     AW          =  ACTIVITY OF WATER
C     B0          =  PITZER PARAMETER, CATION-ANION INTERACTION
C     B1          =  PITZER PARAMETER, CATION-ANION INTERACTION
C     B2          =  PITZER PARAMETER, CATION-ANION INTERACTION
C     BG          =  PITZER PARAMETER, CATION-ANION INTERACTION
C     BGP         =  PITZER PARAMETER, CATION-ANION INTERACTION
C     BM          =  NEUTRAL SPECIE CONCENTRATION
C     BP          =  PITZER PARAM. CATION-ANION INTERACTION
C     CG          =  PITZER PARAMETER, CATION-ANION INTERACTION
C     CI          =  IONIC STRENGTH
C     CM          =  CATION CONCENTRATION
C     CMI         =  IONIC STRENGTH, MOLARITY SCALE
C     CP          =  PITZER PARAMETER, CATION-ANION INTERACTION
C     CS          =  CONC. OF A SPECIES
C     DC          =  DIELECTRIC CONSTANT OF WATER
C     DCA         =  PITZER PARAMETER, NEUTRAL INTERACTION
C     DEN         =  DENSITY OF SOLUTION (G/CM3)
C     DTM         =  PITZER CONST. IN D. H. TERM, (= 1.2)
C     DW          =  DENSITY OF WATER
C     ETHA        =  PITZER PARAMETER FOR UNSYMMETRICAL MIXING OF ANIONS  
C     ETHAP       =  PITZER PARAMETER FOR UNSYMMETRICAL MIXING OF ANIONS
C     ETHC        =  PITZER PARAMETER FOR UNSYMMETRICAL MIXING OF CATIONS  
C     ETHCP       =  PITZER PARAMETER FOR UNSYMMETRICAL MIXING OF CATIONS
C     FG,FGI,FGIP =  FUNCTIONS OF IONIC STRENGTH
C     FP,FPI      =  FUNCTIONS OF IONIC STRENGTH
C     G           =  ACTIVITY COEFFICIENT OF A SPECIES
C     GLNA        =  LN GAMMA FOR ANION 
C     GLNC        =  LN GAMMA FOR CATION 
C     GLNN        =  LN GAMMA FOR NEUTRAL SPECIES
C     IAN         =  INDEX NO. OF AN ANION
C     ICA         =  INDEX NO. OF A CATION
C     INU         =  INDEX NO. OF A NEUTRAL SPECIES
C     NC, NA, NN  =  NO. OF CATIONS, ANIONS, AND NEUTRAL SPECIES
C     NIA         =  INDICATOR OF CALCULATED (0) OR ASSUMED UNIT (1) VALUES
C                     FOR AQU.ACT.COEFS.
C     NUM         =  FLAG INDICATING THAT EFFECTS OF UNSYMMETRICAL MIXING
C                     OF IONS WILL (1) OR WILL NOT (0) BE INCLUDED IN THE
C                     PITZER TREATMENT
C     NZA         =  ARRAY FOR STORING CHARGES ON SPECIES
C     PAA         =  PITZER PARAMETER, ANION-ANION INTERACTION
C     PCAA        =  PITZER PARAMETER, CATION-ANION-ANION INTERACTION
C     PCC         =  PITZER PARAMETER, CATION-CATION INTERACTION)
C     PCCA        =  PITZER PARAMETER, CATION-CATION-ANION INTERACTION
C     PHI         =  THE OSMOTIC COEFFICIENT OF THE SOLN.
C     PNA         =  PITZER PARAMETER, NEUTRAL-ANION INTERACTION
C     PNC         =  PITZER PARAMETER, NEUTRAL-CATION INTERACTION
C     PNN         =  PITZER PARAMETER, NEUTRAL-NEUTRAL INTERACTION
C     P0          =  VAPOR PRESSURE OF WATER
C     RMC         =  RATIO OF MOLALITY TO GIVEN CONCENTRATION
C     TWSA        =  WEIGHT OF SOLUTE (G/L.SLN.,G/KG H2O., OR G/KG SLN.)
C     SV          =  MASSON COEFFICIENT OF AN IONIC SPECIES
C     TM          =  TOTAL MOLALITY OF ALL AQ.SPECIES
C     TN          =  TOTAL EQUIV. CONC. OF ALL CATIONS OR ANIONS
C     TVSA        =  APPARENT VOLUME OF SOLUTE (CM3/L.SLN.,CM3/KG SLV.,
C                       CM3/KG SLN.)
C     TVC0        =  APPARENT VOLUME OF SOLUTE AT INFINITE DILUTION
C                       (CM3/L.SLN.,CM3/KG SLV.,CM3/KG SLN.)
C     V0          =  LIMITING MOLAR VOLUME OF A SPECIES (CM3/MOL)
C     W           =  M.W. OF A SPECIES
C     WA          =  MOLECULAR WEIGHT OF WATER (KG/MOL)
C     ZA          =  ABS. VALUE OF ANION CHARGE
C     ZC          =  ABS. VALUE OF CATION CHARGE
C     ZFA         =  ANION CHARGE FACTOR
C     ZFC         =  CATION CHARGE FACTOR
C
C
      SUBROUTINE AQUAC(P,W,V,V0,SV,CS,AM,
     &BM,CM,B0,B1,B2,CP,ALP,BP,BG,BGP,CG,PCC,ETHC,ETHCP,
     &PAA,ETHA,ETHAP,PCCA,PCAA,DCA,PNC,PNA,PNN,ZA,ZFA,ZC,ZFC,G,
     &ICA,IAN,INU,NZA,ICN,NCAL)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      DIMENSION P(*),W(*),V(*),V0(*),SV(*),CS(*),AM(*),
     &BM(*),CM(*),B0(*),B1(*),B2(*),CP(*),ALP(2,*),
     &BP(*),BG(*),BGP(*),CG(*),PCC(*),ETHC(*),ETHCP(*),
     &PAA(*),ETHA(*),ETHAP(*),PCCA(*),PCAA(*),DCA(*),PNC(*),
     &PNA(*),PNN(*),ZA(*),ZFA(*),ZC(*),ZFC(*),G(*)
      DIMENSION ICA(*),IAN(*),INU(*),NZA(*)
      DIMENSION FPI(2),FGI(2),FGIP(2)
      COMMON/MD/NSO,NCBO,IREF,NPH,NPA,NSA,NIO,NIA,NC,NA,NN,NUP,NUM
      COMMON/DT/TC,NO,NOF,ICNR      
      COMMON/WTR/P0,DW,DC,AP,RLKW      
      COMMON/EQ/DEN,AWG,AW,TWSA,TVSA,VWPA,VWPO,RMBA,RMBO
      DATA TK,DTM,WA/273.15D+0,1.2D+0,0.01801528D+0/
      DATA ZERO,ONE,TWO,THREE,THOU/
     &0.0D+0,1.0D+0,2.0D+0,3.0D+0,1.0D+3/
C  DEFINE FUNCTIONS FOR LOCATING PITZER COEFFS.IN LINEAR ARRAYS
      IX2(I,J,IM)=IM*(J-1)+I
      IX3(I,J,K,IM,JM)=IM*(JM*(K-1)+J-1)+I
C  ASSIGN PITZER PARAMETERS FROM ARRAY P
      IF(NC.GT.0.AND.NA.GT.0) THEN
C      FOR CATION-ANION INTERACTIONS
        K=NSO*(1+NPH+NPA)+2*(NPH+NPA)+NSA+2+NSO
	IJ=0
        DO 15,IC=1,NC
        DO 10 IA=1,NA
        L=IX2(IC,IA,NC)
	IJ=IJ+1
        B0(L)=P(K+1)
        B1(L)=P(K+2)
        B2(L)=P(K+3)
        CP(L)=P(K+4)
        K=K+4
C
C        WRITE(*,1111) IC,IA,L,IJ,B0(L),B1(L),B2(L),CP(L),ALP(1,IJ),
C     &  ALP(2,IJ)
C1111    FORMAT(/' IC IA  L IJ      B0(L)     B1(L)     B2(L)',
C     &  '     CP(L) ALP(1,IJ) ALP(2,IJ)'/4I3,6F10.5)
C        PAUSE
C
10      CONTINUE
15      CONTINUE
      END IF
      IF(NC.GT.1) THEN
C      FOR CATION-CATION INTERACTIONS
        DO 30 IC=1,(NC-1)
        DO 25 JC=(IC+1),NC
        K=K+1
        L=IX2(IC,JC,NC)
        PCC(L)=P(K)
        M=IX2(JC,IC,NC)
        PCC(M)=PCC(L)
        DO 20 IA=1,NA
        K=K+1
        LA=IX3(IC,JC,IA,NC,NC)
        PCCA(LA)=P(K)
        M=IX3(JC,IC,IA,NC,NC)
        PCCA(M)=PCCA(LA)
C
C        WRITE(*,1112) IC,JC,IA,PCC(L),PCCA(M)
C1112    FORMAT(/' IC JC IA    PCC(L)   PCCA(M)'/3I3,2F10.5)
C        PAUSE
C
20      CONTINUE
25      CONTINUE
30      CONTINUE
      END IF
      IF(NA.GT.1) THEN
C      FOR ANION-ANION INTERACTIONS
        DO 50 IA=1,(NA-1)
        DO 45 JA=(IA+1),NA
        K=K+1
        L=IX2(IA,JA,NA)      
        PAA(L)=P(K)
        M=IX2(JA,IA,NA)      
        PAA(M)=PAA(L)
        DO 40 IC=1,NC
        K=K+1
        LC=IX3(IC,IA,JA,NC,NA)
        PCAA(LC)=P(K)
        M=IX3(IC,JA,IA,NC,NA)
        PCAA(M)=PCAA(LC)
C
C        WRITE(*,1113) IC,IA,JA,PAA(L),PCAA(M)
C1113    FORMAT(/' IC IA JA    PAA(L)   PCAA(M)'/3I3,2F10.5)
C        PAUSE
C
40      CONTINUE
45      CONTINUE
50      CONTINUE
      END IF
      IF(NN.GT.0) THEN
C      FOR NEUTRAL-SPECIE INTERACTIONS
        DO 70 IN=1,NN
        K=K+1
        DCA(IN)=P(K)
	IF(NC.GT.0) THEN
          DO 60 IC=1,NC
          K=K+1
          L=IX2(IN,IC,NN) 
          PNC(L)=P(K)
C
C        WRITE(*,1114) IN,IC,L,PNC(L)
C1114    FORMAT(/' IN IC  L    PNC(L)'/3I3,F10.5)
C        PAUSE
C
60        CONTINUE
        END IF
	IF(NA.GT.0) THEN
          DO 65 IA=1,NA
          K=K+1
          L=IX2(IN,IA,NN)
          PNA(L)=P(K)
C
C        WRITE(*,1115) IN,IA,L,PNA(L)
C1115    FORMAT(/' IN IA  L    PNA(L)'/3I3,F10.5)
C        PAUSE
C
65        CONTINUE
        END IF
70      CONTINUE
        DO 80 IN=1,NN
        DO 75 JN=IN,NN
        K=K+1
        L=IX2(IN,JN,NN)
        PNN(L)=P(K)
        M=IX2(JN,IN,NN)
        PNN(M)=PNN(L)
C
C        WRITE(*,1116) IN,JN,L,M,PNN(M)
C1116    FORMAT(/' IN JN  L  M    PNN(M)'/4I3,F10.5)
C        PAUSE
C
75      CONTINUE
80      CONTINUE
      END IF
      T=TC+TK
C  ASSIGN ION CHARGES, VALENCE FACTORS, AND INDICES
C   FOR CATIONIC REACTANTS
      IC=0
      DO 90 I=1,3
      IF(NZA(I).GT.0) THEN
        IC=IC+1
        ZC(IC)=NZA(I)
        ZFC(IC)=DSQRT(ZC(IC)*(ONE+ZC(IC))/TWO)
        IF(I.EQ.1) ICA(IC)=3
        IF(I.EQ.2) ICA(IC)=4
        IF(I.EQ.3) ICA(IC)=7
      END IF
90    CONTINUE
C   FOR ANIONIC REACTANTS
      IA=0
      DO 95 I=4,6
      IF(NZA(I).GT.0) THEN
        IA=IA+1
        ZA(IA)=NZA(I)
        ZFA(IA)=DSQRT(ZA(IA)*(ONE+ZA(IA))/TWO)
        IF(I.EQ.4) IAN(IA)=5
        IF(I.EQ.5) IAN(IA)=6
        IF(I.EQ.6) IAN(IA)=8
      END IF
95    CONTINUE
      IN=0
      IF(NSA.GT.0) THEN
C      FOR AQUEOUS PRODUCT SPECIES
        DO 100 I=1,NSA
        J=8+NSO+I
        K=6+I
        IF(NZA(K).GT.0) THEN
	  IC=IC+1
          ZC(IC)=NZA(K)
          ZFC(IC)=DSQRT(ZC(IC)*(ONE+ZC(IC))/TWO)
          ICA(IC)=J
	END IF
        IF(NZA(K).LT.0) THEN
	  IA=IA+1
          ZA(IA)=-NZA(K)
          ZFA(IA)=DSQRT(ZA(IA)*(ONE+ZA(IA))/TWO)
          IAN(IA)=J
	END IF
        IF(NZA(K).EQ.0) THEN
          IN=IN+1
          INU(IN)=J
	END IF
100     CONTINUE
      END IF
C  VERIFY THE NUMBER OF CATIONS, ANIONS, AND NEUTRALS
      IF(IC.NE.NC.OR.IA.NE.NA.OR.IN.NE.NN) THEN
        WRITE(9,1001)
1001  FORMAT(1X,'CHECK THE NUMBER OF CATIONS. ANIONS, NEUTRALS')
        PAUSE
        STOP
      END IF
      RETURN
C  ENTRY POINT FOR CALCULATING CURRENT ACTIVITY COEFFICIENTS
      ENTRY ENTAAC(P,W,V,V0,SV,CS,AM,
     &BM,CM,B0,B1,B2,CP,ALP,BP,BG,BGP,CG,PCC,ETHC,ETHCP,
     &PAA,ETHA,ETHAP,PCCA,PCAA,DCA,PNC,PNA,PNN,ZA,ZFA,ZC,ZFC,G,
     &ICA,IAN,INU,NZA,ICN,NCAL)
C  ASSIGN CATION, ANION, AND NEUTRAL-SPECIE CONCENTRATIONS;
C   CALCULATE THE IONIC STRENGTH; CALCULATE GRAMS SOLUTE PER LITER
C   SLN., PER KG SLV., OR PER KG SLN.; CALCULATE TOTAL SPECIE,
C   EQUIVALENT, AND IONAL CONCENTRATIONS
      TCI=ZERO
      TWSA=ZERO
      TMC=ZERO
      TNC=ZERO
      IF(NC.GT.0) THEN
        DO 110 I=1,NC
        CM(I)=CS(ICA(I))
        TCI=TCI+CM(I)*ZC(I)**2
        TWSA=TWSA+CM(I)*W(ICA(I))
        TMC=TMC+CM(I)
        TNC=TNC+CM(I)*ZC(I)
110     CONTINUE
      END IF
      IF(NA.GT.0) THEN
        DO 115 I=1,NA
        AM(I)=CS(IAN(I))
        TCI=TCI+AM(I)*ZA(I)**2
        TWSA=TWSA+AM(I)*W(IAN(I))
        TMC=TMC+AM(I)
115     CONTINUE
      END IF
      IF(NN.GT.0) THEN
        DO 120 I=1,NN
        BM(I)=CS(INU(I))
        TWSA=TWSA+BM(I)*W(INU(I))
        TMC=TMC+BM(I)
120     CONTINUE
      END IF
      CMI=TCI/TWO
      IF(ICN.EQ.1) THEN
C      SUM THE APPARENT MOLAR VOLUME OF EACH SPECIES (AT THE MOLAR 
C       IONIC STRENGTH OF THE SOLUTION) TIMES ITS CONCENTRATION
        SCMI=DSQRT(CMI)
        TVSA=ZERO
	IF(NC.GT.0) THEN
          DO 130 I=1,NC
          IF(ICA(I).LT.9) K=ICA(I)-2
          IF(ICA(I).GT.8) K=ICA(I)-2-NSO
          TVSA=TVSA+(V0(K)+SV(K)*SCMI/ZFC(I))*CM(I)
130       CONTINUE
        END IF
	IF(NA.GT.0) THEN
          DO 135 I=1,NA
          IF(IAN(I).LT.9) K=IAN(I)-2
          IF(IAN(I).GT.8) K=IAN(I)-2-NSO
          TVSA=TVSA+(V0(K)+SV(K)*SCMI/ZFA(I))*AM(I)
135       CONTINUE
        END IF
        IF(NN.GT.0) THEN
          DO 140 I=1,NN
          K=INU(I)-2-NSO
          TVSA=TVSA+V0(K)*BM(I)
140       CONTINUE
        END IF
C      CALCULATE THE DENSITY OF THE SOLUTION
        DEN=DW*(ONE-TVSA/THOU)+TWSA/THOU
C      CALCULATE THE RATIO OF MOLALITY TO THE GIVEN CONC.
        RMC=ONE/(DEN-TWSA/THOU)
      END IF
      IF(ICN.EQ.2) RMC=ONE
      IF(ICN.EQ.3) RMC=ONE/(ONE-TWSA/THOU)
      IF(NIA.EQ.1) THEN
C      UNIT SOLUTE ACTIVITY COEFFICIENTS ARE ASSUMED
        IF(NC.GT.0) THEN
          DO 155 I=1,NC
          G(ICA(I))=ONE
155       CONTINUE
        END IF
        IF(NA.GT.0) THEN
          DO 160 I=1,NA
          G(IAN(I))=ONE
160       CONTINUE
        END IF
        IF(NN.GT.0) THEN
          DO 165 I=1,NN
          G(INU(I))=ONE
165       CONTINUE
        END IF
        IF(AWG.GT.ZERO) THEN
C        THE WATER ACTIVITY HAS BEEN SUPPLIED BY THE USER
          AW=AWG
	  RETURN
        END IF
C      CALCULATE WATER ACTIVITY
        IF(ICN.EQ.1) THEN
          TVC0=ZERO
	  IF(NC.GT.0) THEN
            DO 175 I=1,NC
            IF(ICA(I).LT.9) K=ICA(I)-2
            IF(ICA(I).GT.8) K=ICA(I)-2-NSO
            TVC0=TVC0+V0(K)*CM(I)
175         CONTINUE
          END IF
	  IF(NA.GT.0) THEN
            DO 180 I=1,NA
            IF(IAN(I).LT.9) K=IAN(I)-2
            IF(IAN(I).GT.8) K=IAN(I)-2-NSO
            TVC0=TVC0+V0(K)*AM(I)
180         CONTINUE
          END IF
          IF(NN.GT.0) THEN
            DO 185 I=1,NN
            K=INU(I)-2-NSO
            TVC0=TVC0+V0(K)*BM(I)
185         CONTINUE
          END IF
	  AW=DEXP((THOU*WA*TMC/(DW*TVC0))*DLOG(ONE-TVC0/THOU))
	END IF
        IF(ICN.EQ.2) THEN
	  AW=DEXP(-TMC*WA)
	END IF
        IF(ICN.EQ.3) THEN
	  AW=DEXP((THOU*WA*TMC/TWSA)*DLOG(ONE-TWSA/THOU))
	END IF
        RETURN	  
      END IF
C  ACTIVITY COEFFICIENTS ARE TO BE CALCULATED BY THE PITZER TREATMENT
C   CONVERT CONCENTRATIONS TO MOLALITIES
      IF(NC.GT.0) THEN
        DO 200 I=1,NC
        CM(I)=CM(I)*RMC
200     CONTINUE
      END IF
        IF(NA.GT.0) THEN
        DO 205 I=1,NA
        AM(I)=AM(I)*RMC
205     CONTINUE
      END IF
      IF(NN.GT.0) THEN
        DO 210 I=1,NN
        BM(I)=BM(I)*RMC
210     CONTINUE
      END IF
      TM=TMC*RMC
      TN=TNC*RMC
      CI=CMI*RMC
C  CALCULATE NEEDED FUNCTIONS OF THE IONIC STRENGTH
      SCI=DSQRT(CI)
      F1I=ONE+DTM*SCI
      FP=-AP*SCI/F1I
      FG=FP-AP*TWO/DTM*DLOG(F1I)
      IF(NC.GT.0.AND.NA.GT.0) THEN
C      CALCULATE B(I,J) AND C(I,J) PARAMETERS
        IJ=0
        DO 235 I=1,NC
        DO 230 J=1,NA
        L=IX2(I,J,NC)
	IJ=IJ+1
        DO 225 K=1,2
        FPI(K)=DEXP(-ALP(K,IJ)*SCI)
        IF(CI.GT.ZERO.AND.ALP(K,IJ).GT.ZERO) THEN 
          FGI(K)=(ONE-(ONE+ALP(K,IJ)*SCI)*FPI(K))/(ALP(K,IJ)**2*CI)
          FGIP(K)=(-ONE+(ONE+ALP(K,IJ)*SCI+0.5*ALP(K,IJ)**2*CI)*
     &    FPI(K))/(ALP(K,IJ)*CI)**2
        END IF
        IF(CI.EQ.ZERO.OR.ALP(K,IJ).EQ.ZERO) THEN 
          FGI(K)=0.5
C         IF ALP(2,IJ) IS 0 THEN FGIP(2) SHOULD BE SET AT 0
C         IF CI IS 0 THEN FGIP GOES TO INFINITY, BUT ITS PRODUCT WITH
C          CONC**2 GOES TO ZERO.  AS CI --> 0, FGIP --> -ALP/(2*SCI)
          FGIP(K)=ZERO
        END IF
225     CONTINUE
        BP(L)=B0(L)+B1(L)*FPI(1)+B2(L)*FPI(2)
        BG(L)=B0(L)+TWO*B1(L)*FGI(1)+TWO*B2(L)*FGI(2)
        BGP(L)=TWO*B1(L)*FGIP(1)+TWO*B2(L)*FGIP(2)        
        CG(L)=CP(L)/(TWO*DSQRT(ZC(I)*ZA(J)))
230     CONTINUE
235     CONTINUE
      END IF
      IF(NUM.EQ.1) THEN 
        IF(NC.GT.1) THEN
C        EVALUATE FUNCTIONS NEEDED FOR UNSYMM. MIXING OF CATIONS
          DO 238 I=1,NC-1
          DO 237 J=I+1,NC
          L=IX2(I,J,NC)
          M=IX2(J,I,NC)
          ETHC(L)=ZERO
          ETHCP(L)=ZERO
          IF(ZC(I).NE.ZC(J).AND.CI.GT.ZERO)
     &    CALL ETHETA(CI,AP,ZC(I),ZC(J),ETHC(L),ETHCP(L))
          ETHC(M)=ETHC(L)
          ETHCP(M)=ETHCP(L)
237       CONTINUE
238       CONTINUE
        END IF
        IF(NA.GT.1) THEN
C        EVALUATE FUNCTIONS NEEDED FOR UNSYMM. MIXING OF ANIONS
          DO 242 I=1,NA-1
          DO 241 J=I+1,NA
          L=IX2(I,J,NA)
          M=IX2(J,I,NA)
          ETHA(L)=ZERO
          ETHAP(L)=ZERO
          IF(ZA(I).NE.ZA(J).AND.CI.GT.ZERO)
     &    CALL ETHETA(CI,AP,ZA(I),ZA(J),ETHA(L),ETHAP(L))
          ETHA(M)=ETHA(L)
          ETHAP(M)=ETHAP(L)
241       CONTINUE
242       CONTINUE
        END IF
      END IF
C  CALCULATE THE OSMOTIC COEFFICIENT AND WATER ACTIVITY
      IF(TM.GT.ZERO) CALL PHIMIX(ZC,ZA,CM,AM,BM,TM,TN,
     &CI,FP,BP,CP,PCC,ETHC,ETHCP,PAA,ETHA,ETHAP,PCCA,PCAA,DCA,PNC,
     &PNA,PNN,PHI)
      IF(TM.EQ.ZERO) PHI=ONE
      AW=DEXP(-TM*WA*PHI)      
C  CALCULATE ACTIVITY COEFFICIENTS AND CONVERT TO USER CONC. SCALE
      IF(NC.GT.0) THEN
        DO 240 I=1,NC
        CALL CGAMA(ZC,ZA,CM,AM,BM,TN,FG,
     &  BG,BGP,CG,PCC,ETHC,ETHCP,ETHAP,PCCA,PCAA,DCA,PNC,I,GLNC)
        K=ICA(I)
        IF(ICN.EQ.1) G(K)=RMC*DW*DEXP(GLNC)
        IF(ICN.EQ.2) G(K)=DEXP(GLNC)
        IF(ICN.EQ.3) G(K)=RMC*DEXP(GLNC)
240     CONTINUE
      END IF
      IF(NA.GT.0) THEN
        DO 245 I=1,NA
        CALL AGAMA(ZC,ZA,CM,AM,BM,TN,FG,
     &  BG,BGP,CG,PAA,ETHA,ETHAP,ETHCP,PCCA,PCAA,DCA,PNA,I,GLNA)
        K=IAN(I)
        IF(ICN.EQ.1) G(K)=RMC*DW*DEXP(GLNA)
        IF(ICN.EQ.2) G(K)=DEXP(GLNA)
        IF(ICN.EQ.3) G(K)=RMC*DEXP(GLNA)
245     CONTINUE
      END IF
      IF(NN.GT.0) THEN
        DO 250 I=1,NN
        CALL NGAMA(ZC,ZA,CM,AM,BM,TN,DCA,
     &  PNC,PNA,PNN,I,GLNN)
        K=INU(I)
        IF(ICN.EQ.1) G(K)=RMC*DW*DEXP(GLNN)
        IF(ICN.EQ.2) G(K)=DEXP(GLNN)
        IF(ICN.EQ.3) G(K)=RMC*DEXP(GLNN)
250     CONTINUE
      END IF
      RETURN
      END
C
C
C
      SUBROUTINE WATER(TC,P0,DW,DC,RLKW)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      DATA ONE,ZERO,DIFF/1.0D+0,0.0D+0,1.0D-10/
      CALL PSAT(TC,P0)
      DW=ONE
      DO 20 I=1,100
      CALL WAPVT(TC,DW,PX,BETA)
      XS=DW
      DW=DW/DEXP(BETA*(PX-P0))
      IF(DABS((DW-XS)/XS).LT.DIFF) GO TO 21
20    CONTINUE
      WRITE(9,100) 
100   FORMAT(1X,'CALCULATED WATER DENSITY NOT CORRECT')
      PAUSE
      STOP
21    CALL WDBP(TC,DW,P0,DC)
      CALL KWMF(TC,DW,RLKW)
      RETURN
      END
C
C
C
C  SUBROUTINE TO CALCULATE THE SATURATION PRESSURE OF PURE WATER
C   BASED ON KEENAN ET AL (1969), EQ 17
C
      SUBROUTINE PSAT(TC,P)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      DIMENSION F(8)
      DATA TK,TCRT,PCRT/273.15D+0,374.136D+0,220.88D+0/
      DATA F/-741.9242D+0,-29.72100D+0,-11.55286D+0,-8.685635D-1,
     &1.094098D-1,4.39993D-1,2.520658D-1,5.218684D-2/
      DATA ONE,ZERO,CN1,CN2,CN3,CN4/1.0D+0,0.0D+0,1.0D+3,0.65D+0,
     &0.01D+0,1.0D-5/
      T=TC+TK
      TR=CN1/T
      S=ZERO
      DO 10 I=1,8
      S=S+F(I)*(CN2-CN3*TC)**(I-1)
10    CONTINUE
      P=PCRT*DEXP(CN4*TR*(TCRT-TC)*S)
      RETURN
      END
C
C
C
C  SUBROUTINE TO CALCULATE THE PROPERTIES OF PURE WATER
C   FROM THE STEAM TABLES OF KEENAN ET AL (1969)
C
C    X = DENSITY (G/CM3), P = PRESSURE (BAR), TC = THERMO. TEMP. (C),
C
      SUBROUTINE WAPVT(TC,X,P,BETA)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      DIMENSION C(8),A(10,7),TI(8)
      DATA TK,TCR,E,RP/273.15D+0,1.544912D+0,4.8D+0,4.6151D+0/
      DATA C/1.857065D+3,3.22912D+3,-4.19465D+2,3.66649D+1,-2.05516D+1
     &,4.85233D+0,4.6D+1,-1.011249D+3/
      DATA A/2.9492937D+1,-1.3213917D+2,2.7464632D+2,-3.6093828D+2
     &,3.4218431D+2,-2.4450042D+2,1.5518535D+2,5.9728487D+0
     &,-4.1030848D+2,-4.1605860D+2
     &,-5.1985860D+0,7.7779182D+0,-3.3301902D+1,-1.6254622D+1
     &,-1.7731074D+2,1.2748742D+2,1.3746153D+2,1.5597836D+2
     &,3.3731180D+2,-2.0988866D+2
     &,6.8335354D+0,-2.6149751D+1,6.5326396D+1,-2.6181978D+1
     &,0.0D+0,0.0D+0,0.0D+0,0.0D+0,-1.3746618D+2,-7.3396848D+2
     &,-1.564104D-1,-7.2546108D-1,-9.2734289D+0,4.3125840D+0
     &,0.0D+0,0.0D+0,0.0D+0,0.0D+0,6.7874983D+0,1.0401717D+1
     &,-6.3972405D+0,2.6409282D+1,-4.7740374D+1,5.6323130D+1
     &,0.0D+0,0.0D+0,0.0D+0,0.0D+0,1.3687317D+2,6.4581880D+2
     &,-3.9661401D+0,1.5453061D+1,-2.9142470D+1,2.9568796D+1
     &,0.0D+0,0.0D+0,0.0D+0,0.0D+0,7.9847970D+1,3.9917570D+2
     &,-6.9048554D-1,2.7407416D+0,-5.1028070D+0,3.9636085D+0
     &,0.0D+0,0.0D+0,0.0D+0,0.0D+0,1.3041253D+1,7.1531353D+1/
      DATA ONE,ZERO,TEN,CN1,CN2,RMW,CN4/1.0D+0,0.0D+0,10.0D+0,
     &1.0D+3,2.0D+0,18.0153D+0,0.634D+0/
      RC=RP/TEN
      T=TC+TK
      TR=CN1/T
      TF=TR-TCR
      FK=DEXP(-E*X)
      SC=ZERO
      DO 10 I=1,6
      SC=SC+C(I)/TR**(I-1)
10    CONTINUE
      AF0=SC+C(7)*DLOG(T)+C(8)*DLOG(T/TR)
      Q=ZERO
      DQDX=ZERO
      DQDX2=ZERO
      DO 20 J=1,7
      IF(J.GT.1) GO TO 21
      TP=TCR
      XP=CN4
      GO TO 22
21    TP=2.5D+0
      XP=ONE
22    XT=X-XP
      TT=TR-TP
      TI(1)=ONE
      DO 30 I=2,8
      TI(I)=TI(I-1)*XT
30    CONTINUE
      SI=A(1,J)+A(2,J)*TI(2)+A(3,J)*TI(3)+A(4,J)*TI(4)+A(5,J)*TI(5)
     &+A(6,J)*TI(6)+A(7,J)*TI(7)+A(8,J)*TI(8)
      DIDX=A(2,J)+2.0*A(3,J)*TI(2)+3.0*A(4,J)*TI(3)+4.0*A(5,J)*TI(4)+
     &5.0*A(6,J)*TI(5)+6.0*A(7,J)*TI(6)+7.0*A(8,J)*TI(7)
      DIDX2=2.0*A(3,J)+6.0*A(4,J)*TI(2)+12.0*A(5,J)*TI(3)+20.0*A(6,J)*
     &TI(4)+30.0*A(7,J)*TI(5)+42.0*A(8,J)*TI(6)
      SK=FK*(A(9,J)+A(10,J)*X)
      DKDX=FK*A(10,J)-E*SK
      DKDX2=-E*FK*A(10,J)-E*DKDX
      TJ=TT**(J-2)
      Q=Q+TJ*(SI+SK)
      DQDX=DQDX+TJ*(DIDX+DKDX)
      DQDX2=DQDX2+TJ*(DIDX2+DKDX2)
20    CONTINUE
      Q=TF*Q
      DQDX=TF*DQDX
      DQDX2=TF*DQDX2
      AF=AF0+RC*T*(DLOG(X)+X*Q)
      P=X*RP*T*(ONE+X*Q+X**2*DQDX)
      DPDX=P/X+X*RP*T*(Q+X*DQDX+CN2*X*DQDX+X**2*DQDX2)
      BETA=ONE/(X*DPDX)
      GG=AF+(P/X)*(RC/RP)
      G=GG*RMW
      RETURN
      END
C
C
C
C  SUBROUTINE TO CALCULATE THE DIELECTRIC CONSTANT OF WATER
C   BASED ON EQ (1) OF BRADLEY AND PITZER (1979)
C
      SUBROUTINE WDBP(TC,DW,P0,DC)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      DIMENSION U(9)
      DATA U/3.4279D+2,-5.0866D-3,9.4690D-7,-2.0525D+0,3.1159D+3,
     &-1.8289D+2,-8.0325D+3,4.2142D+6,2.1417D+0/
      DATA ONE,ZERO,CN1/1.0D+0,0.0D+0,1.0D+3/
      TK=273.15D+0
      T=TC+TK
      DCO=U(1)*DEXP(U(2)*T+U(3)*T**2)
      C=U(4)+U(5)/(U(6)+T)
      B=U(7)+U(8)/T+U(9)*T
      P=P0
      DC=DCO+C*DLOG((B+P)/(B+CN1))
      DT1=U(1)*(U(2)+2.0*U(3)*T)*DEXP(U(2)*T+U(3)*T**2)
      DT2=(U(5)/(U(6)+T)**2)*DLOG((B+P)/(B+CN1))
      DT3=(C*(CN1-P)*(U(9)-U(8)/T**2))/((B+P)*(B+CN1))
      DCDT=DT1-DT2+DT3
      RETURN
      END
C
C
C  SUBROUTINE TO CALCULATE ETHETA AND D(ETHETA)/D(CI) FOR
C   UNSYMMETRICAL MIXING OF LIKE-CHARGED IONS
C
C
      SUBROUTINE ETHETA(CI,AP,Z1,Z2,ETH,ETHP)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 J11,J12,J22,J11P,J12P,J22P
      DATA HALF,FOUR,SIX,EIGHT/0.5D+00,4.0D+00,
     &6.0D+00,8.0D+00/
      FA=SIX*AP*DSQRT(CI)
      X12=Z1*Z2*FA
      X11=Z1*Z1*FA
      X22=Z2*Z2*FA
      CALL UNSYM(X11,J11,J11P)
C      WRITE(13,200) X11,J11,J11P
      CALL UNSYM(X12,J12,J12P)
C      WRITE(13,210) X12,J12,J12P
      CALL UNSYM(X22,J22,J22P)
C      WRITE(13,220) X22,J22,J22P
      ETH=(Z1*Z2/(FOUR*CI))*(J12-HALF*J11-HALF*J22)
      ETHP=-(ETH/CI)+(Z1*Z2/(EIGHT*CI**2))*
     &(X12*J12P-HALF*X11*J11P-HALF*X22*J22P)
      ETHPI=ETHP*CI
C      WRITE(13,230) CI,ETH,ETHPI
C200   FORMAT(1X,'XMM, JMM, JMMP = ',3F14.7)
C210   FORMAT(1X,'XMN, JMN, JMNP = ',3F14.7)
C220   FORMAT(1X,'XNN, JNN, JNNP = ',3F14.7/)
C230   FORMAT(1X,'I, ETHETA, CI*DETHDI = ', 3F9.4/)
      RETURN
      END
C
C
C
C  SUBROUTINE TO CALCULATE KW OF WATER FROM THE MARSHALL-FRANK EQUATION
C
C         TC = TEMP / C
C         DW = DENSITY OF WATER AT SPECIFIED T AND P=I ATM
C       RLKW = LOG10(KW)
C
C
      SUBROUTINE KWMF(TC,DW,RLKW)
      IMPLICIT REAL*8(A-H,O-Z)
      DATA ONE,THOU/1.0D+00,1.0D+03/
      DATA A,B,C,D/
     &-4.35284D+00,-2978.479D+00,1.33968D+05,-3.03364D+07/
      DATA E,F,G/13.957D+00,-1262.3D+00,8.5641D+05/
      DATA T0,R,RLN/273.15D+00,8.314510D-03,2.30258509D+00/
      T=TC+T0
      F1=A+B/T+C/T**2+D/T**3
      F2=E+F/T+G/T**2
      RLKW = F1+F2*DLOG10(DW)
      RETURN
      END
C
C
C
C   SUBROUTINE TO CALCULATE J AND J' FOR UMSYMMETRICAL MIXING OF IONS
C
C     TREATMENT IS FROM PITZER V, EQ. 47, P. 263
C     (J. SOLN. CHEM., 4, 249-265, 1975)
C
C          X = Z(I)*Z(J)*E**2/(DKT)
C          J = F(X)
C       DJDX = DJ/DX
C
C
      SUBROUTINE  UNSYM(X,J,DJDX)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 J,DJDX
      DATA C1,C2,C3,C4,ONE,FOUR/
     &4.581D+00,0.7237D+00,0.0120D+00,0.528D+00,1.0D+00,4.0D+00/
      XF=C1*X**(-C2)
      YF=-C3*X**C4
      EF=DEXP(YF)
      J=X/(FOUR+XF*EF)
      DXFDX=-C2*XF/X
      DYFDX=C4*YF/X
      DEFDX=DYFDX*EF
      DFDX=DXFDX*EF+XF*DEFDX
      DJDX=J/X-J**2*DFDX/X
      RETURN
      END
C
C
C
C
C
C
C  SUBROUTINE TO CALCULATE THE OSMOTIC COEFFICIENT OF A MIXTURE
C   FROM THE TREATMENT OF PITZER
C
      SUBROUTINE PHIMIX(ZC,ZA,CM,AM,BM,TM,TN,CI,FP,BP,CP,PCC,
     &ETHC,ETHCP,PAA,ETHA,ETHAP,PCCA,PCAA,DCA,PNC,PNA,PNN,PHI)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      COMMON/MD/NSO,NCBO,IREF,NPH,NPA,NSA,NIO,NIA,NC,NA,NN,NUP,NUM
      COMMON/WTR/P0,DW,DC,AP,RLKW      
      DIMENSION ZC(*),ZA(*),CM(*),AM(*),BM(*),BP(*),CP(*),
     &PCC(*),ETHC(*),ETHCP(*),PAA(*),ETHA(*),ETHAP(*),PCCA(*),
     &PCAA(*),DCA(*),PNC(*),PNA(*),PNN(*)
      DATA ZERO,ONE,TWO/0.0D+00,1.0D+00,2.0D+00/
C  DEFINE FUNCTIONS FOR LOCATING PITZER COEFFS.IN LINEAR ARRAYS
      IX2(I,J,IM)=IM*(J-1)+I
      IX3(I,J,K,IM,JM)=IM*(JM*(K-1)+J-1)+I
C
      TP1=TWO*CI*FP
      TP2=ZERO
      DO 10 I=1,NC
      DO 15 J=1,NA
      L=IX2(I,J,NC)
      TP2=TP2+TWO*CM(I)*AM(J)*(BP(L)+(TN/DSQRT(ZC(I)*ZA(J)))*
     &CP(L))
15    CONTINUE
10    CONTINUE
      TP3=ZERO
      IF(NC.LT.2) GO TO 35
      DO 20 I=1,NC
      DO 25 J=1,NC
      IF(I.EQ.J) GO TO 25
      SMA=ZERO
      DO 30 K=1,NA
      L=IX3(I,J,K,NC,NC)
      SMA=SMA+AM(K)*PCCA(L)
30    CONTINUE
      L=IX2(I,J,NC)
      PCCP=PCC(L)
      IF(NUM.EQ.1.AND.ZC(I).NE.ZC(J)) PCCP=PCCP+ETHC(L)+CI*ETHCP(L)
      TP3=TP3+CM(I)*CM(J)*(PCCP+SMA)
25    CONTINUE
20    CONTINUE
35    TP4=ZERO
      IF(NA.LT.2) GO TO 60
      DO 40 I=1,NA
      DO 45 J=1,NA
      IF(I.EQ.J) GO TO 45
      SMC=ZERO
      DO 50 K=1,NC
      L=IX3(K,I,J,NC,NA)
      SMC=SMC+CM(K)*PCAA(L)
50    CONTINUE
      L=IX2(I,J,NA)
      PAAP=PAA(L)
      IF(NUM.EQ.1.AND.ZA(I).NE.ZA(J)) PAAP=PAAP+ETHA(L)+CI*ETHAP(L)
      TP4=TP4+AM(I)*AM(J)*(PAAP+SMC)
45    CONTINUE
40    CONTINUE
60    TP5=ZERO
      IF(NN.EQ.0) GOTO 90
      DO 85 I=1,NN
      SMC=ZERO
      DO 70 J=1,NC
      M=IX2(I,J,NN)
      SMC=SMC+ZC(J)*CM(J)*PNC(M)
70    CONTINUE
      SMA=ZERO
      DO 75 K=1,NA
      M=IX2(I,K,NN)
      SMA=SMA+ZA(K)*AM(K)*PNA(M)
75    CONTINUE
      DO 80 L=1,NN
      SMB=ZERO
      M=IX2(I,L,NN)
      SMB=SMB+BM(L)*PNN(M)
80    CONTINUE
      TP5=TP5+BM(I)*(TN*DCA(I)+SMC+SMA+SMB)
85    CONTINUE
90    PHI=ONE+(TP1+TP2+TP3+TP4+TP5)/TM
      RETURN
      END
C
C
C
C  SUBROUTINE TO CALCULATE THE ACTIVITY COEFFICIENT OF A CATION
C   FROM THE TREATMENT OF PITZER
C
      SUBROUTINE CGAMA(ZC,ZA,CM,AM,BM,TN,FG,BG,BGP,CG,PCC,
     &ETHC,ETHCP,ETHAP,PCCA,PCAA,DCA,PNC,I,GLNC)
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      COMMON/MD/NSO,NCBO,IREF,NPH,NPA,NSA,NIO,NIA,NC,NA,NN,NUP,NUM
      COMMON/WTR/P0,DW,DC,AP,RLKW      
      DIMENSION ZC(*),ZA(*),CM(*),AM(*),BM(*),BG(*),BGP(*),
     &CG(*),PCC(*),ETHC(*),ETHCP(*),ETHAP(*),PCCA(*),PCAA(*),
     &DCA(*),PNC(*)
      DATA ZERO,HALF,ONE,TWO/0.0D+0,0.5D+0,1.0D+0,2.0D+0/
C  DEFINE FUNCTIONS FOR LOCATING PITZER COEFFS.IN LINEAR ARRAYS
      IX2(ID,J,IM)=IM*(J-1)+ID
      IX3(ID,J,K,IM,JM)=IM*(JM*(K-1)+J-1)+ID
C
      TG1=ZC(I)**2*FG
      TG2=ZERO
      DO 5 J=1,NA
      L=IX2(I,J,NC)
      TG2=TG2+TWO*AM(J)*(BG(L)+TN*CG(L))
5     CONTINUE    
      TG3=ZERO
      DO 10 J=1,NC
      DO 15 K=1,NA
      L=IX2(J,K,NC)
      TG3=TG3+CM(J)*AM(K)*(ZC(I)**2*BGP(L)+ZC(I)*CG(L))
15    CONTINUE
10    CONTINUE
      TG4=ZERO
      DO 20 J=1,NC
      IF(J.EQ.I) GO TO 20
      SMA=ZERO
      DO 25 K=1,NA
      L=IX3(I,J,K,NC,NC)
      SMA=SMA+AM(K)*PCCA(L)
25    CONTINUE
      L=IX2(I,J,NC)      
      PCCP=PCC(L)
      IF(NUM.EQ.1.AND.ZC(I).NE.ZC(J)) PCCP=PCCP+ETHC(L)
      TG4=TG4+CM(J)*(TWO*PCCP+SMA)
20    CONTINUE
      TG5=ZERO