c http://Evgenii.Rudnyi.Ru/
c
      PROGRAM DKCl
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      EXTERNAL SS
      PARAMETER (IXJAC = 600)
      DIMENSION X(4), PARM(4), F(600), XJAC(IXJAC,4), XJTJ(10),
     *          WORK(1230), XtXinv(4)
      CHARACTER Name*25, CurTime*8
      LOGICAL EQ
      COMMON /Dat/ NEXT, NT, NEXTR, NTR, NEXEF, NEF, NEXK5, NK5,
     *    NAME(40), NP(40), NFP(40), T(600), P(600), P1(600), EQ(40)
      COMMON /Work/ AinvT(600), HS1(600), HS2(600), R
      DIMENSION AvrInvT(40), PT(40)
      R = 8.31441
      CALL SetKCl1(.True.)

      DO i = 1, NTr
        AinvT(i) = 1e3/T(i)
        CALL HSCOR (T(i), DH1, DS1, DH2, DS2)
        HS1(i) = (-DH1/T(i) + DS1)/R
        HS2(i) = EXP((-DH2/T(i) + DS2)/R)
      END DO
      DO i = 1, NexTr
        AvrInvT(i) = 0.
        DO j = NFP(i), NFP(i) + NP(i) - 1
          AvrInvT(i) = AvrInvT(i) + AinvT(j)
        END DO
        AvrInvT(i) = AvrInvT(i)/NP(i)
        PT(i) = 0.
        DO j = NFP(i), NFP(i) + NP(i) - 1
          PT(i) = PT(i) + (AinvT(j) - AvrInvT(i))**2
        END DO
      END DO

      OPEN (3, NAME='kcl.lst')
      CALL UGETIO (3, 5, 3)
      N = 4

C   JANAF
      X(1) = 206.466
      X(2) = 132.139
      X(3) = 23.645
      X(4) = 16.252
      CALL SS(X, NTr, N, F)
      SSQ = 0.
      DO i = 1, NTr
        SSQ = SSQ + F(i)**2
      END DO
      DO i = 1, NexTr
        IF (Eq(i)) SSQ = SSQ + 0.03**2*(NP(i)-2)
      END DO
      WRITE (3, '(4F11.3)') X
      WRITE (3, *) 'Sum of squares for JANAF ', SSQ

C   Gurvich
      X(1) = 207.150
      X(2) = 133.030
      X(3) = 26.257
      X(4) = 14.483
      CALL SS(X, NTr, N, F)
      SSQ = 0.
      DO i = 1, NTr
        SSQ = SSQ + F(i)**2
      END DO
      DO i = 1, NexTr
        IF (Eq(i)) SSQ = SSQ + 0.03**2*(NP(i)-2)
      END DO
      WRITE (3, '(4F11.3)') X
      WRITE (3, *) 'Sum of squares for Gurvich ', SSQ
      WRITE (3, *)

C  ML(III_1)
      X(1) = 206.775
      X(2) = 133.030
      X(3) = 37.834
      X(4) = 24.913
      CALL SS(X, NTr, N, F)
      SSQ = 0.
      DO i = 1, NTr
        SSQ = SSQ + F(i)**2
      END DO
      DO i = 1, NexTr
        IF (Eq(i)) SSQ = SSQ + 0.03**2*(NP(i)-2)
      END DO
      WRITE (3, *) 'Sum of squares for ML(III-1)', SSQ
      StDev = SQRT(SSQ/(NTr - N))
      WRITE (3, 1) X, SSQ, NTr, StDev
      WRITE (3, *)
      WRITE (3, 5)
      WRITE (3, *) 'Total pressure'
      DO i = 1, NexTr
        S = 0.
        Sa = 0.
        Sb = 0.
        DO j = NFP(i), NFP(i) + NP(i) - 1
          S = S + F(j)**2
          Sa = Sa + F(j)
          Sb = Sb + F(j)*(AinvT(j) - AvrInvT(i))
        END DO
        Sa = Sa/NP(i)
        IF (PT(i) .GT. 0.) Sb = Sb/PT(i)
        Se1 = S - Sa**2*NP(i) - Sb**2*PT(i)
        IF (NP(i) .GT. 2) THEN
          Ste1 = Se1/(NP(i)-2)
        ELSE
          Ste1 = 0.
        END IF
        WRITE (3, 4) Name(i), i,
     *    SQRT(S/NP(i)), SQRT(Ste1), Sa, Sb, Sb*SQRT(PT(i)/NP(i)),
     *    Sa-(0.8/5.5*0.1), (Sb*SQRT(PT(i)/NP(i)))-(0.2/5.*0.05)
        IF (i .EQ. NexT) THEN
          WRITE (3, *)
          WRITE (3, *) 'Knudsen effusion'
        ELSE IF (i .EQ. NexEf) THEN
          WRITE (3, *)
          WRITE (3, *) 'Transpiration'
        END IF
      END DO
4     FORMAT (A16, I3, 7F7.3)
5     FORMAT ('                  i  sqrt(S/NP) sqrt(Ste1)',
     * 'Sa Sb Sb1')

c      X(1) = 207.003
c      X(2) = 132.535
c      X(3) = 34.201
c      X(4) = 23.928
      X(1) = 207.228
      X(2) = 133.061
      X(3) = 36.290
      X(4) = 24.252
      NSIG = 4
      EPS = 0.
      DELTA = 0.
      MAXFN = 100
      IOPT = 1
      CALL TIME(CurTime)
      PRINT*, CurTime
      WRITE (3, *) CurTime
      PRINT*, 'Starting ZXSSQ'

      CALL ZXSSQ(SS,NTr,N,NSIG,EPS,DELTA,MAXFN,IOPT,
     *           PARM,X,SSQ,F,XJAC,IXJAC,XJTJ,WORK,INFER,IER)
      WRITE (3, *) 'The norm of the gradient ', WORK(1)
      WRITE (3, *) 'The number of function evaluation ', WORK(2)
      WRITE (3, *) 'The number of significant digits in parameters '
     *                   , WORK(3)
      WRITE (3, *) 'The MARQUARDT scaling parameter', WORK(4)
      WRITE (3, *) 'The number of iteration ', WORK(5)
      WRITE (3, *) 'INFER = ', INFER, '  IER = ', IER
      DO i = 1, NexTr
        IF (Eq(i)) SSQ = SSQ + 0.03**2*(NP(i)-2)
      END DO
      StDev = SQRT(SSQ/(NTr - N))
      WRITE (3, 1) X, SSQ, NTr, StDev
1     FORMAT (/,1X, '    DelH1     DelS1     DelH2     DelS2',/,
     *  1X, 4F11.3,//, 'Sum of squares',F10.4, '  Number of points',
     *  I5,/, 'Standard deviation [sqrt(SS/(N-4))]', F10.3,/)

      IJOB = 1
      CALL LINV3P (XJTJ, PARM, IJOB, N, IER)
      WRITE (3, *) 'LINV3P - IER ', IER
      WRITE (3, *) '   (XJ''*XJ)**(-1)'
      DO i = 1, N
        WRITE (3, '(+1P, 4E10.2)')(XJTJ(j), j = i*(i-1)/2+1, i*(i+1)/2)
      END DO
      WRITE (3, *)
      WRITE (3, *) '   Parameter  sqrt(Cii)  st. dev.'
      DO i = 1, N
        XtXinv(i) = SQRT(XJTJ(i*(i+1)/2))
        WRITE (3, '(4F11.2)') X(i), XtXinv(i), XtXinv(i)*StDev
      END DO
      WRITE (3, *)
      WRITE (3, *) 'Correlation matrix'
      WRITE (3, '(F10.3)') 1
      DO i = 2, N
        DO j = 1, i-1
          Loc = i*(i-1)/2+j
          XJTJ(Loc) = XJTJ(Loc)/SQRT(XJTJ(i*(i+1)/2)*XJTJ(j*(j+1)/2))
        END DO
        WRITE (3, '(4F10.3)') (XJTJ(j), j = i*(i-1)/2+1, i*(i+1)/2-1), 1
      END DO

      CALL TIME(CurTime)
      PRINT*, CurTime
      WRITE (3, *)
      WRITE (3, *) CurTime
      PRINT*, 'Data used in minimization'
      WRITE (3, *) 'Data used in minimization'
      WRITE (3, 3)
      WRITE (3, *) 'Total pressure'
      SteAvr = 0.
      Nu = 0
      DO i = 1, NexTr
        S = 0.
        Se = 0.
        Sa = 0.
        Sb = 0.
        A = 0.
        B = 0.
        DO j = NFP(i), NFP(i) + NP(i) - 1
          A = A + P(j)
          B = B + P(j)*(AinvT(j) - AvrInvT(i))
          S = S + F(j)**2
          Sa = Sa + F(j)
          Sb = Sb + F(j)*(AinvT(j) - AvrInvT(i))
        END DO
        A = A/NP(i)
        Sa = Sa/NP(i)
        IF (PT(i) .GT. 0.) Sb = Sb/PT(i)
        IF (PT(i) .GT. 0.) B = B/PT(i)
        DO j = NFP(i), NFP(i) + NP(i) - 1
          Se = Se + (P(j) - A - B*(AinvT(j) - AvrInvT(i)))**2
        END DO
        Se1 = S - Sa**2*NP(i) - Sb**2*PT(i)
        SteAvr = SteAvr + Se1
        IF (.NOT. Eq(i)) Nu = Nu + NP(i) - 2
        IF (NP(i) .GT. 2) THEN
          Ste = Se/(NP(i)-2)
          Ste1 = Se1/(NP(i)-2)
        ELSE
          Ste = 0.
          Ste1 = 0.
        END IF
        WRITE (3, 2) Name(i), i,
     *             SQRT(S/NP(i)), SQRT(Ste), SQRT(Ste1), Sa, Sb,
     *             Sb*SQRT(PT(i)/NP(i)), A - B*AvrInvT(i), B
        IF (i .EQ. NexT) THEN
          WRITE (3, *)
          WRITE (3, *) 'Knudsen effusion'
        ELSE IF (i .EQ. NexEf) THEN
          WRITE (3, *)
          WRITE (3, *) 'Transpiration'
        END IF
      END DO
      WRITE (3, *)
      WRITE (3, *) 'SteAvr ' , SQRT(SteAvr/Nu), '  Nu ', Nu
2     FORMAT (A16, I3, 6F7.3, F6.2, F7.2)
3     FORMAT ('                  i  sqrt(S/NP) sqrt(Ste) sqrt(Ste1)',
     * 'Sa Sb Sb1  A    B')

      CLOSE (1)
      END

      SUBROUTINE SS(X, M, N, F)
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      DIMENSION X(N), F(M)
      CHARACTER Name*25
      LOGICAL EQ
      COMMON /Dat/ NEXT, NT, NEXTR, NTR, NEXEF, NEF, NEXK5, NK5,
     *    NAME(40), NP(40), NFP(40), T(600), P(600), P1(600), EQ(40)
      COMMON /Work/ AinvT(600), HS1(600), HS2(600), R
      DO i = 1, NT
        PM = (-X(1)*AinvT(i) + X(2))/R + HS1(i)
        PDM = EXP((-X(3)*AinvT(i) + X(4))/R)*HS2(i)
        F(i) = PM + LOG(1. + PDM) - P(i)
      END DO
      DO i = NT + 1, NEf
        PM = (-X(1)*AinvT(i) + X(2))/R + HS1(i)
        PDM = EXP((-X(3)*AinvT(i) + X(4))/R)*HS2(i)
        F(i) = PM + LOG(1. + SQRT(2.d0)*PDM) - P(i)
      END DO
      DO i = NEf + 1, NTr
        PM = (-X(1)*AinvT(i) + X(2))/R + HS1(i)
        PDM = EXP((-X(3)*AinvT(i) + X(4))/R)*HS2(i)
        F(i) = PM + LOG(1. + 2.*PDM) -
     *      LOG(1. + EXP(PM)*PDM/P1(i))  - P(i)
      END DO
      RETURN
      END