c http://Evgenii.Rudnyi.Ru/
c
      function i_trim(b)
      character b*8
      i_trim = 8
      do while (b(i_trim:i_trim) .eq. ' ' .and. i_trim .ne. 0)
       i_trim = i_trim - 1
      end do
      return
      end

      FUNCTION Locij(i, j)
      IF (i .GE. j) THEN
        Locij = j + (i*i - i)/2
      ELSE
        Locij = i + (j*j - j)/2
      END IF
      RETURN
      END

      PROGRAM EQCONST
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      PARAMETER (ZERO = 0d0, ONE = 1d0, R = 8.31441d0)
      PARAMETER (MaxR = 30, MaxS = 30)
      DOUBLE PRECISION S, EKplus, EKprod, Hcalc, Scalc, Hr, Sr, N
      INTEGER Ms, MHa, MHb, MSa, MSb, Nr, i, j, nwork, k, Ni, jj, j1,
     * jj1, Ijob, Ier, is, if, NCov, Nu
      DIMENSION Hexp(MaxR), Sexp(MaxR), P(MaxR), N(MaxR), T(MaxR)
      DIMENSION X(MaxR,MaxS)
      DIMENSION Hf(MaxS), Sf(MaxS), AMass(MaxS), AIsotope(MaxS),
     *   sHf(MaxS), sSf(MaxS), HT(MaxS), Sp(MaxS), Sf2(MaxS)
      CHARACTER Name(MaxS)*8, Line*80, ChaWork(20)*8, NowTime*8
      LOGICAL MassCor, ExpData, Ion(MaxS)
      DIMENSION Work(20), Temp(50), EK(50), Si(50), Ni(50)
      DIMENSION A(2*MaxS,2*MaxR), C((2*MaxR)*(2*MaxR+1)/2)
      DIMENSION A1(2*MaxS,2*MaxS)
      DIMENSION D(2*MaxS,2*MaxS), Cov(2*MaxS,2*MaxS)
      DIMENSION Covkl(40), kCov(40), lCov(40)

      CALL GETARG(1, Line)
      is = 1
      DO WHILE (Line(is:is) .EQ. ' ' .AND. is .LT. 70)
        is = is + 1
      END DO
      i = 0
      if = is
      DO WHILE (Line(if:if) .NE. ' ')
       IF (Line(if:if) .EQ. '.') i = if
       if = if + 1
      END DO
      IF (i .EQ. 0) THEN
        i = if
        if = if + 4
        Line(i:if-1) = '.dat'
      END IF
      OPEN (1, NAME=Line(is:if-1))
      OPEN (2, NAME=Line(is:i)//'lst')
      CALL TIME(NowTime)
      PRINT*, NowTime
      WRITE (2, *) NowTime
      READ (1, 100) Line
      DO WHILE (Line(1:1) .EQ. ';' .OR. LEN_TRIM(Line) .EQ. 0)
        READ (1, 100) Line
      END DO
      j = 0
      DO WHILE (LEN_TRIM(Line) .NE. 0)
        j = j + 1
        READ (Line, *) Name(j), Ion(j), AMass(j), AIsotope(j),
     *      Hf(j), HT(j), sHf(j), Sf(j), sSf(j)
        READ (1, 100) Line
      END DO
      Ms = j
      j = 1
      MHa = j
      DO WHILE (Hf(j) .EQ. ZERO)
        j = j + 1
        MHa = j
      END DO
      j = 1
      MSa = j
      DO WHILE (Sf(j) .EQ. ZERO)
        j = j + 1
        MSa = j
      END DO
      MHb = MHa - 1
      MSb = MSa - 1

      WRITE (2,'(I4, '' substances'', /, I4, '' substances with unknown
     *enthalpies'', /, I4, '' substances with unknown entropies'')')
     * Ms, MHb, MSb
      WRITE (2, *) 'name    ion    M       i     DelfH0    HT-H0',
     *  ' s(H)   ST   s(S)'
      DO j = 1, Ms
        IF (j .EQ. MHa) WRITE (2, *) 'Substances with known enthalpies'
        IF (j .EQ. MSa) WRITE (2, *) 'Substances with known entropies'
        WRITE (2, '(A, L3, F7.1, F8.4, F10.1, F7.1, F5.1, F7.1, F5.1)')
     *   Name(j), Ion(j), AMass(j), AIsotope(j), Hf(j), HT(j), sHf(j),
     *   Sf(j), sSf(j)
      END DO
      WRITE (2, *)
      CALL TIME(NowTime)
      PRINT*, NowTime
      WRITE (2, *) NowTime

C  read covariations H and H
      READ (1, 100) Line
      DO WHILE (Line(1:1) .EQ. ';' .OR. LEN_TRIM(Line) .EQ. 0)
        READ (1, 100) Line
      END DO
      i = 0
      IF (Line(1:1) .EQ. '*') GO TO 3
      WRITE (2, *)
      WRITE (2, *) 'Covariations H and H'
      DO WHILE (LEN_TRIM(Line) .NE. 0)
        i = i + 1
        READ (Line, *) Covkl(i), ChaWork(1), ChaWork(2)
        DO j = MHa, Ms
          IF (Name(j) .EQ. ChaWork(1)) THEN
            kCov(i) = j - MHa + 1
            GO TO 1
          END IF
        END DO
        STOP 'There is no such substance'
1       CONTINUE
        DO j = MHa, Ms
          IF (Name(j) .EQ. ChaWork(2)) THEN
            lCov(i) = j - MHa + 1
            GO TO 2
          END IF
        END DO
        STOP 'There is no such substance'
2       CONTINUE
        WRITE (2, '(F8.2, 2X, A, I2, 1X, A, I2)')
     *    Covkl(i), ChaWork(1), kCov(i), ChaWork(2), lCov(i)
        READ (1, 100) Line
      END DO

C  read covariations H and S
3     READ (1, 100) Line
      DO WHILE (Line(1:1) .EQ. ';' .OR. LEN_TRIM(Line) .EQ. 0)
        READ (1, 100) Line
      END DO
      IF (Line(1:1) .EQ. '*') GO TO 6
      WRITE (2, *)
      WRITE (2, *) 'Covariations H and S'
      DO WHILE (LEN_TRIM(Line) .NE. 0)
        i = i + 1
        READ (Line, *) Covkl(i), ChaWork(1), ChaWork(2)
        DO j = MHa, Ms
          IF (Name(j) .EQ. ChaWork(1)) THEN
            kCov(i) = j - MHa + 1
            GO TO 4
          END IF
        END DO
        STOP 'There is no such substance'
4       CONTINUE
        DO j = MSa, Ms
          IF (Name(j) .EQ. ChaWork(2)) THEN
            lCov(i) = j - MSa + 1 + Ms - MHa + 1
            GO TO 5
          END IF
        END DO
        STOP 'There is no such substance'
5       CONTINUE
        WRITE (2, '(F8.2, 2X, A, I2, 1X, A, I2)')
     *    Covkl(i), ChaWork(1), kCov(i), ChaWork(2), lCov(i)
        READ (1, 100) Line
      END DO

C  read covariations S and S
6     READ (1, 100) Line
      DO WHILE (Line(1:1) .EQ. ';' .OR. LEN_TRIM(Line) .EQ. 0)
        READ (1, 100) Line
      END DO
      IF (Line(1:1) .EQ. '*') GO TO 16
      WRITE (2, *)
      WRITE (2, *) 'Covariations S and S'
      DO WHILE (LEN_TRIM(Line) .NE. 0)
        i = i + 1
        READ (Line, *) Covkl(i), ChaWork(1), ChaWork(2)
        DO j = MHa, Ms
          IF (Name(j) .EQ. ChaWork(1)) THEN
            kCov(i) = j - MSa + 1 + Ms - MHa + 1
            GO TO 14
          END IF
        END DO
        STOP 'There is no such substance'
14      CONTINUE
        DO j = MSa, Ms
          IF (Name(j) .EQ. ChaWork(2)) THEN
            lCov(i) = j - MSa + 1 + Ms - MHa + 1
            GO TO 15
          END IF
        END DO
        STOP 'There is no such substance'
15      CONTINUE
        WRITE (2, '(F8.2, 2X, A, I2, 1X, A, I2)')
     *    Covkl(i), ChaWork(1), kCov(i), ChaWork(2), lCov(i)
        READ (1, 100) Line
      END DO
16    NCov = i
      WRITE (2, *)

C  read reactions
      i = 1
      DO
        READ (1, 100, END = 10) Line
        DO WHILE (Line(1:1) .EQ. ';' .OR. LEN_TRIM(Line) .EQ. 0)
          READ (1, 100, END = 10) Line
        END DO
        nwork = 0
        is = 1
        DO WHILE (is .LT. LEN_TRIM(Line))
          DO WHILE (Line(is:is) .EQ. ' ')
            is = is + 1
          END DO
          if = is
          DO WHILE (Line(if:if) .NE. ' ')
            if = if + 1
          END DO
        nwork = nwork + 1
        READ (Line(is:if-1), *) ChaWork(nwork)
        is = if
        END DO
        READ (1, *) (Work(k), k = 1, nwork)
        READ (1, *) Sum, MassCor, ExpData
        DO j = 1, Ms
          X(i,j) = ZERO
        END DO
        DO k = 1, nwork
          DO j = 1, Ms
            IF (Name(j) .EQ. ChaWork(k)) THEN
              X(i,j) = Work(k)
              GO TO 8
            END IF
          END DO
          STOP 'There is no such substance'
8         CONTINUE
        END DO
        S = 0.
        DO j = 1, Ms
          IF (Ion(j)) S = S + X(i,j)**2
        END DO
        EKprod = ONE
        IF (Sum .EQ. ZERO) THEN
          Sum = SQRT(S)
          EKprod = ONE/Sum
          WRITE (2, *)
     *       'Equilibrium constants will be divided by SQRT(Sum)'
        ELSE IF (Sum .NE. S) THEN
          WRITE (2, *) 'Reaction is not normalised, squared norm is', S
          Sum = SQRT(Sum)
        ELSE
          Sum = SQRT(Sum)
        END IF
        DO j = 1, Ms
          X(i,j) = X(i,j)/Sum
        END DO
        WRITE (2, '(''('',I2,'')'',SP,20(F5.0, 1X, A))')
     *     i, (Work(k), ChaWork(k)(1:i_trim(ChaWork(k))), k = 1, nwork)
        WRITE (2, '(''Sum is '', F8.0)') Sum**2
        EKplus = ZERO
        IF (MassCor) THEN
          DO j = 1, Ms
            IF(Ion(j))
     *         EKplus = EKplus + X(i,j)*LOG(AMass(j)/AIsotope(j))
          END DO
          WRITE (2, '(''Mass correction will be '', F7.3)') EKplus
        END IF
        WRITE (2, *)
        WRITE (2, *) ' Tavr      b      a   DelH            DelS
     *     N      P     Sp    Sf2'
        WRITE (2, *)
        PRINT*, 'Reaction N ', i

C  read and process equilibrium constants
        READ (1, 100) Line
        DO WHILE (Line(1:1) .EQ. ';' .OR. LEN_TRIM(Line) .EQ. 0)
          READ (1, 100) Line
        END DO
        IF (ExpData) THEN
          if = 0
          DO WHILE (Line(1:1) .NE. '*')
            is = if + 1
            Hexp(i) = ZERO
            Sexp(i) = ZERO
            P(i) = ZERO
            T(i) = ZERO
            N(i) = ZERO
            Sp(i) = ZERO
            Sf2(i) = ZERO
            DO WHILE (LEN_TRIM(Line) .NE. 0)
              if = if + 1
              READ (Line, *) Temp(if), EK(if), Si(if), Ni(if)
              EK(if) = EK(if)*EKprod + EKplus
              Si(if) = Si(if)*EKprod
              WRITE (2, '(F6.0, F8.2, F6.2, I4)')
     *              Temp(if), EK(if), Si(if), Ni(if)
              Temp(if) = 1000d0/Temp(if)
              T(i) = T(i) + Temp(if)*Ni(if)
              Sexp(i) = Sexp(i) + EK(if)*Ni(if)
              N(i) = N(i) + Ni(if)
              IF (Ni(if) .GT. 1)
     *            Sp(i) = Sp(i) + Si(if)**2*(Ni(if)-1)
              READ (1, 100) Line
            END DO
            T(i) = T(i)/N(i)
            Sexp(i) = Sexp(i)/N(i)
            IF (N(I) .GT. if-is+1) Sp(i) = Sp(i)/(N(i)-if+is-1)
            DO k = is, if
              P(i) = P(i) + (Temp(k) - T(i))**2*Ni(k)
              Hexp(i) = Hexp(i) + EK(k)*(Temp(k) - T(i))*Ni(k)
            END DO
            Hexp(i) = Hexp(i)/P(i)
            DO k = is, if
              Sf2(i) = Sf2(i) +
     *          (EK(k) - Hexp(i)*(Temp(k) - T(i)) - Sexp(i))**2*Ni(k)
            END DO
            IF (if .GT. is+1) Sf2(i) = Sf2(i)/(if-is-1)
            WRITE (2, *)
            WRITE (2,'(F6.0, 2F7.2,
     *          2(F7.1,'' +/-'',F5.1),F5.1,F8.3,2F6.2)')
     *         1000d0/T(i), Hexp(i), Sexp(i),
     *         -Hexp(i)*R, SQRT(Sf2(i)/P(i))*R,
     * (Sexp(i)-Hexp(i)*T(i))*R, SQRT(Sf2(i)*(T(i)**2/P(i)+One/N(i)))*R,
     *         N(i), SQRT(P(i)), SQRT(Sp(i)), SQRT(Sf2(i))
            WRITE (2, *)
            DO WHILE (Line(1:1) .EQ. ';' .OR. LEN_TRIM(Line) .EQ. 0)
              READ (1, 100, END = 10) Line
            END DO
            i = i + 1
            DO j = 1, Ms
              X(i,j) = X(i-1,j)
            END DO
          END DO
          WRITE (2, *) 'All experiments'
          Hexp(i) = ZERO
          Sexp(i) = ZERO
          P(i) = ZERO
          T(i) = ZERO
          N(i) = ZERO
          Sp(i) = ZERO
          Sf2(i) = ZERO
          DO k = 1, if
            T(i) = T(i) + Temp(k)*Ni(k)
            Sexp(i) = Sexp(i) + EK(k)*Ni(k)
            N(i) = N(i) + Ni(k)
            IF (Ni(k) .GT. 1)
     *          Sp(i) = Sp(i) + Si(k)**2*(Ni(k)-1)
          END DO
          T(i) = T(i)/N(i)
          Sexp(i) = Sexp(i)/N(i)
          IF (N(i) .GT. if) Sp(i) = Sp(i)/(N(i)-if)
          DO k = 1, if
            P(i) = P(i) + (Temp(k) - T(i))**2*Ni(k)
            Hexp(i) = Hexp(i) + EK(k)*(Temp(k) - T(i))*Ni(k)
          END DO
          Hexp(i) = Hexp(i)/P(i)
          DO k = 1, if
            Sf2(i) = Sf2(i) +
     *          (EK(k) - Hexp(i)*(Temp(k) - T(i)) - Sexp(i))**2*Ni(k)
          END DO
          IF (if .GT. 2) Sf2(i) = Sf2(i)/(if-2)
            WRITE (2,'(F6.0, 2F7.2,
     *          2(F7.1,'' +/-'',F5.1),F5.1,F8.3,2F6.2)')
     *         1000d0/T(i), Hexp(i), Sexp(i),
     *         -Hexp(i)*R, SQRT(Sf2(i)/P(i))*R,
     * (Sexp(i)-Hexp(i)*T(i))*R, SQRT(Sf2(i)*(T(i)**2/P(i)+One/N(i)))*R,
     *         N(i), SQRT(P(i)), SQRT(Sp(i)), SQRT(Sf2(i))
        ELSE
          READ (Line, *) T(i), Hexp(i), Sexp(i),
     *                       N(i), P(i), Sp(i), Sf2(i)
          T(i) = 1000d0/T(i)
          Hexp(i) = -Hexp(i)/R
          Sexp(i) = Sexp(i)/R + Hexp(i)*T(i)
          P(i) = P(i)**2
          Sp(i) = Sp(i)**2
          Sf2(i) = Sf2(i)**2
          IF (P(i) .NE. Zero) THEN
            Hr = SQRT(Sf2(i)/P(i))*R
          ELSE
            Hr = 0.
          END IF
          IF (N(i) .NE. Zero) THEN
            Sr = SQRT(Sf2(i)/(P(i)*T(i)**2+One/N(i)))*R
          ELSE
            Sr = SQRT(Sf2(i)/P(i)*T(i)**2)*R
          END IF
          WRITE (2,'(F6.0, 2(F7.1,'' +/-'',F5.1),F5.1,F8.3,2F6.2)')
     *         1000d0/T(i), -Hexp(i)*R, Hr,
     *         (Sexp(i)-Hexp(i)*T(i))*R, Sr,
     *         N(i), SQRT(P(i)), SQRT(Sp(i)), SQRT(Sf2(i))
          i = i + 1
        END IF
        WRITE (2, *)
      END DO
10    Nr = i - 1
      WRITE (2, *) 'Overall there are ', Nr, ' reactions'
      WRITE (2, *)
      CALL TIME(NowTime)
      PRINT*, NowTime
      WRITE (2, *) NowTime
      CLOSE (1)

C  modifying Hexp and Sexp
      DO i = 1, Nr
        Hexp(i) = Hexp(i)*R
        DO j = MHa, Ms
          Hexp(i) = Hexp(i) + X(i,j)*(Hf(j)+HT(j))
        END DO
      END DO
      DO i = 1, Nr
        Sexp(i) = Sexp(i)*R
        DO j = MHa, Ms
          Sexp(i) = Sexp(i) + T(i)*X(i,j)*(Hf(j)+HT(j))
        END DO
        DO j = MSa, Ms
          Sexp(i) = Sexp(i) - X(i,j)*Sf(j)
        END DO
      END DO

C  computing (X'D^(-1)X)^(-1)
      DO j = 1, MHb
        jj = (j*j - j)/2
        DO j1 = 1, j
          jj1 = j1 + jj
          C(jj1) = Zero
          DO i = 1, Nr
            C(jj1) = C(jj1) + X(i,j)*(P(i)+T(i)*N(i)*T(i))*X(i,j1)
          END DO
        END DO
      END DO
      DO j = MHb + 1, MHb + MSb
        jj = (j*j - j)/2
        DO j1 = 1, MHb
          jj1 = j1 + jj
          C(jj1) = Zero
          DO i = 1, Nr
            C(jj1) = C(jj1) - X(i,j-MHb)*N(i)*T(i)*X(i,j1)
          END DO
        END DO
        DO j1 = MHb + 1, j
          jj1 = j1 + jj
          C(jj1) = Zero
          DO i = 1, Nr
            C(jj1) = C(jj1) + X(i,j-MHb)*N(i)*X(i,j1-MHb)
          END DO
        END DO
      END DO
      Ijob = 1
c      CALL LINV3P (C, A, Ijob, MHb + MSb, Ier)
      call dpptrf('U', MHb + MSb, C, Ier)
      IF (Ier. GT. 0) STOP 'DPPTRF can''t work'
      call dpptri('U', MHb + MSb, C, Ier)
      IF (Ier. GT. 0) STOP 'DPPTRI can''t work'

      WRITE (2, *)
      CALL TIME(NowTime)
      PRINT*, 'Inverting is over ', NowTime
      WRITE (2, *) 'Inverting is over ', NowTime

C  computing a transformation matrix (X'D^(-1)X)^(-1) X'D^(-1)
      DO i = 1, Nr
        DO j = 1, MHb + MSb
          A(j,i) = Zero
          DO j1 = 1, MHb
            A(j,i) = A(j,i) - C(Locij(j,j1))*X(i,j1)*P(i)
          END DO
        END DO
      END DO
      DO i = Nr + 1, 2*Nr
        DO j = 1, MHb
          A(j,i) = Zero
          DO j1 = 1, MHb
            A(j,i) = A(j,i) - C(Locij(j,j1))*X(i-Nr,j1)*T(i-Nr)*N(i-Nr)
          END DO
          DO j1 = MHb + 1, MHb + MSb
            A(j,i) = A(j,i) + C(Locij(j,j1))*X(i-Nr,j1-MHb)*N(i-Nr)
          END DO
        END DO
      END DO
      DO i = Nr + 1, 2*Nr
        DO j = MHb + 1, MHb + MSb
          A(j,i) = Zero
          DO j1 = 1, MHb
            A(j,i) = A(j,i) - C(Locij(j,j1))*X(i-Nr,j1)*T(i-Nr)*N(i-Nr)
          END DO
          DO j1 = MHb + 1, MHb + MSb
            A(j,i) = A(j,i) + C(Locij(j,j1))*X(i-Nr,j1-MHb)*N(i-Nr)
          END DO
        END DO
      END DO

C  computing a solution
      DO j = 1, MHb
        Hf(j) = Zero
        DO i = 1, Nr
          Hf(j) = Hf(j) + A(j,i)*Hexp(i)
        END DO
        DO i = 1, Nr
          Hf(j) = Hf(j) + A(j,i+Nr)*Sexp(i)
        END DO
      END DO
      DO j = 1, MSb
        Sf(j) = Zero
        DO i = 1, Nr
          Sf(j) = Sf(j) + A(j+MHb,i)*Hexp(i)
        END DO
        DO i = 1, Nr
          Sf(j) = Sf(j) + A(j+MHb,i+Nr)*Sexp(i)
        END DO
      END DO

      WRITE (2, *)
      CALL TIME(NowTime)
      PRINT*, 'Solution is found ', NowTime
      WRITE (2, *) 'Solution is found ', NowTime

C  computing a sum of squares and printing reactions
      WRITE (2, *)
      WRITE (2, *) ' i  DelHcor  DelHest  dif  DelScor DelSest  dif
     * DelH  DelS   S1'
      S = Zero
      DO i = 1, Nr
        Hcalc = Zero
        Scalc = Zero
        DO j = 1, MHb
          Hcalc = Hcalc - X(i,j)*Hf(j)
          Scalc = Scalc - T(i)*X(i,j)*Hf(j)
        END DO
        DO j = 1, MSb
          Scalc = Scalc + X(i,j)*Sf(j)
        END DO
        Hr = Zero
        Sr = Zero
        DO j = 1, MHb
          Hr = Hr + X(i,j)*Hf(j)
        END DO
        DO j = MHa, Ms
          Hr = Hr + X(i,j)*(Hf(j) + HT(j))
        END DO
        DO j = 1, Ms
          Sr = Sr + X(i,j)*Sf(j)
        END DO
        S1 = (Hexp(i) - Hcalc)**2*P(i) + (Sexp(i) - Scalc)**2*N(i)
        WRITE (2, '(I3, 2(2F8.1, F6.1), 2F8.1, F6.2)')
     *  i, Hexp(i), Hcalc, Hexp(i)-Hcalc, Sexp(i), Scalc,
     *  Sexp(i)-Scalc, Hr, Sr, SQRT(S1)/R
        S = S + S1
      END DO
      Nu = 2*Nr - MHb - MSb
      DO i = 1, Nr
        IF (N(i) .EQ. Zero) Nu = Nu - 1
        IF (P(i) .EQ. Zero) Nu = Nu - 1
      END DO
      IF (Nu .GT. 0) S = S/Nu

      WRITE (2, *)
      WRITE (2, '(A, F6.2)') 'Standart error of fitting is ', SQRT(S)/R
      WRITE (2, '(A, I3)') 'Degrees of freedom are ', Nu
      WRITE (2, *)
      S = S*4

C  computing A1 - to estimate variance and covariance matrices
      DO j1 = 1, Ms-MHa+1
        DO j = 1, MHb+MSb
          A1(j,j1) = Zero
          DO i = 1, Nr
            A1(j,j1) = A1(j,j1) - A(j,i)*X(i, j1+MHa-1)
          END DO
          DO i = Nr+1, 2*Nr
            A1(j,j1) = A1(j,j1) - A(j,i)*X(i-Nr, j1+MHa-1)*T(i-Nr)
          END DO
        END DO
      END DO
      DO j1 = Ms-MHa+2, Ms-MHa+1+Ms-MSa+1
        DO j = 1, MHb+MSb
          A1(j,j1) = Zero
          DO i = Nr + 1, 2*Nr
            A1(j,j1) = A1(j,j1) + A(j,i)*X(i-Nr, j1-(Ms-MHa+2)+MSa)
          END DO
        END DO
      END DO

C  computing covariance matrix
      DO j1 = 1, Ms-MHa+1+Ms-MSa+1
        DO j = 1, MHb+MSb
          Cov(j, j1) = Zero
        END DO
      END DO
      DO j1 = 1, Ms-MHa+1
        DO j = 1, MHb+MSb
          Cov(j, j1) = Cov(j, j1) + A1(j, j1)*sHf(j1+MHa-1)**2
        END DO
      END DO
      DO j1 = Ms-MHa+2, Ms-MHa+1+Ms-MSa+1
        DO j = 1, MHb+MSb
          Cov(j, j1) = Cov(j, j1) + A1(j, j1)*sSf(j1-(Ms-MHa+2)+MSa)**2
        END DO
      END DO
      DO j1 = 1, NCov
        DO j = 1, MHb+MSb
          Cov(j,lCov(j1)) = Cov(j,lCov(j1)) + A1(j,kCov(j1))*Covkl(j1)
        END DO
        DO j = 1, MHb+MSb
          Cov(j,kCov(j1)) = Cov(j,kCov(j1)) + A1(j,lCov(j1))*Covkl(j1)
        END DO
      END DO

C  computing variance matrix
      DO j1 = 1, MHb+MSb
        DO j = 1, MHb+MSb
          D(j, j1) = Zero
          DO i = 1, Ms-MHa+1+Ms-MSa+1
            D(j, j1) = D(j, j1) + Cov(j, i)*A1(j1, i)
          END DO
        END DO
      END DO

C  printing
      WRITE (2, *)
     * '         DelfHT   DelfH0  sr    sold  ssum'
      DO j = 1, MHb
        sHf(j) = SQRT(S*C(Locij(j,j)) + D(j,j))
        WRITE (2, '(A, 2F8.1, 3F6.1)')
     *     Name(j), Hf(j), Hf(j) - HT(j), SQRT(S*C(Locij(j,j))),
     *     SQRT(D(j,j)), sHf(j)
      END DO
      WRITE (2, *)
      WRITE (2, *) '            S      sr   sold  ssum'
      DO j = 1, MSb
        sSf(j) = SQRT(S*C(Locij(j+MHb,j+MHb)) + D(j+MHb,j+MHb))
        WRITE (2, '(A, F8.1, 3F6.1)')
     *     Name(j), Sf(j), SQRT(S*C(Locij(j+MHb,j+MHb))),
     *     SQRT(D(j+MHb,j+MHb)), sSf(j)
      END DO
      DO j1 = 1, MHb+Msb
        DO j = 1, MHb+Msb
          D(j,j1) = D(j,j1)+S*C(Locij(j,j1))
        END DO
      END DO
      WRITE (2, *)
      WRITE (2, '(8X, 100A6)')
     *   (Name(i), i = 1, MHb), (Name(i), i = 1, MSb)
      DO j = 1, MHb
        WRITE (2, '(A8, 100F6.0)')
     *       Name(j), (D(j,j1), j1 = 1, MHb+MSb)
      END DO
      DO j = MHb+1, MHb+MSb
        WRITE (2, '(A8, 100F6.0)')
     *      Name(j-MHb), (D(j,j1), j1 = 1, MHb+MSb)
      END DO
      WRITE (2, *)
      DO j = 1, Ms-MHa+1
        WRITE (2, '(A8, 100F6.0)')
     *      Name(j+MHa-1), (-Cov(j1,j), j1 = 1, MHb+MSb)
      END DO
      WRITE (2, *)
      DO j = Ms-MHa+2, Ms-MHa+1+Ms-MSa+1
        WRITE (2, '(A8, 100F6.0)')
     *      Name(j-(MS-MHa+2)+MSa), (-Cov(j1,j), j1 = 1, MHb+MSb)
      END DO

      WRITE (2, *)
      CALL TIME(NowTime)
      PRINT*, 'All done  ', NowTime
      WRITE (2, *) 'All done  ', NowTime
      CLOSE (2)

100   FORMAT(A)
      END