PROGRAM MAINPRG
        PARAMETER (MVAR1=8,NTOHOP=400, NNHOM=3)
        INTEGER IDX_YT(12), IDK_YT(12)       ! 12 MONTH IN YEAR
        INTEGER IDX_NT(12), IDK_NT(12)       ! 12 MONTH IN YEAR
        
      REAL PTR(NTOHOP,NNHOM,MVAR1+1), ACCURATE(NTOHOP) ! LUU CAC PHUONG TRINH VA DO CHINH XAC
        CHARACTER HSPTR(NTOHOP,NNHOM,MVAR1+1)*10, VARNAME(MVAR1)*10
        
      INTEGER SHIFT, IPT
        CHARACTER *20 FNAME_NT(8), FNAME_YT
        CHARACTER FOUT*70
      character *12 fsl,fra,FTEMP

        NAMELIST  /FN_IN_OUT/ FSL, FTEMP, FRA 
        1          /CONTROL/   SHIFT, IDX_YT, IDK_YT, IDX_NT, IDK_NT
        2          /FNAMENT/   FNAME_NT
        3          /FNAMEYT/   FNAME_YT, FOUT

        COMMON /STORE_KQ/ PTR, ACCURATE, IPT
        COMMON /FNAME/   fsl,fra,FTEMP
        COMMON /CTRL/    SHIFT,IDX_YT,IDK_YT,IDX_NT,IDK_NT
        COMMON /FNAM_NT/ FNAME_NT
        COMMON /FNAM_YT/ FNAME_YT, FOUT

        DATA VARNAME /
        &       '    NINO12','    NINO_3','    NINO34','    NINO_4',
        &       '       SOI','     SST_1','     SST_2','     SST_3' /

        PTR=-999.0
        HSPTR='      ----'

      OPEN (21,FILE='INPUT_DISR.TXT',STATUS='OLD')
        READ (21,FN_IN_OUT)
        READ (21,CONTROL)
      READ (21,FNAMENT  )
      READ (21,FNAMEYT  )

        CLOSE (21)
        
      open(3,File=FRA,status='unknown')

      CALL TACHSL    !...TACH SO LIEU THANH CAC LOP
      CALL CHON_TINH

        DO I=1,IPT
           DO J=1,NNHOM
              DO K=1,MVAR1+1
                 IF (PTR(I,J,K).NE.-999.0) THEN
                    WRITE(HSPTR(I,J,K),'(F10.4)') PTR(I,J,K)
                 ENDIF
              ENDDO
           ENDDO
        ENDDO

        DCXMAX=0.0
        IPTMAX=0

        OPEN (55,FILE=FOUT,STATUS='UNKNOWN')
        
        WRITE (55,'("     ","   HESO_Co",10A10)') (VARNAME(I),I=1,MVAR1)
      DO I=1,IPT
           WRITE(55,'("PTRINH ",I3)') I
           WRITE(55,'("DCXAC: ",F10.1)') ACCURATE(I)
           DO J=1,NNHOM
!             WRITE (55,'("HAM",I2,10F10.4)') J,(PTR(I,J,K),K=1,MVAR1+1)
              WRITE (55,'("HAM",I2,10A10)') J,(HSPTR(I,J,K),K=1,MVAR1+1)
           ENDDO
           IF (ACCURATE(I).GT.DCXMAX) THEN
              DCXMAX=ACCURATE(I)
              IPTMAX=I
           ENDIF
        ENDDO
        WRITE(55,*) '======================================'
        WRITE(55,'(" P.TRINH THU ",I3," - DO C/XAC MAX = ",F10.4)') 
        &              IPTMAX, DCXMAX


      END ! MAIN PROGRAM
!####################################################

      SUBROUTINE TACHSL
        PARAMETER (NMAX=200, MMAX=12, KMAX=3,MVAR=8, N=40)
        REAL X(NMAX,0:MMAX), XX(NMAX*MMAX,MVAR), YY(NMAX*MMAX)
        REAL XG(NMAX*MMAX,MVAR,KMAX) ! MANG CHUA SO LIEU CAC NHOM
                                     ! (KMAX NHOM)
        INTEGER IDX_YT(12), IDK_YT(12)       ! 12 MONTH IN YEAR
        INTEGER IDX_NT(12), IDK_NT(12)       ! 12 MONTH IN YEAR
        INTEGER SHIFT
        INTEGER NN,YEAR1,YEAR2
        CHARACTER *20 FNAME_NT(MVAR), FNAME_YT
        COMMON /CTRL/    SHIFT,IDX_YT,IDK_YT,IDX_NT,IDK_NT
        COMMON /FNAM_NT/ FNAME_NT
        COMMON /FNAM_YT/ FNAME_YT

C...MAKE CODE
        CALL MAKECODE(FNAME_YT)

C...READ PREDICTANT-DATA 

        OPEN (41,FILE='DRAFT.TXT',STATUS='OLD')
!       OPEN (41,FILE='TEMP\CODE_B1T.DAT',STATUS='OLD')
!       DO I=1,15
        READ(41,*)
        READ(41,*) NN, YEAR1, YEAR2
        READ(41,*)
!      ENDDO

        M=MMAX
        DO I=1,N
           READ(41,*) (X(I,J),J=0,M)
        ENDDO
        CLOSE(1)

        IF (SHIFT.EQ.1) THEN ! (1)
      NN=0
      DO K=1,12
         IF (IDK_YT(K).EQ.0) THEN
                  DO I=1,N-1
                 NN=NN+1
                 YY(NN)=X(I,IDX_YT(K))
              ENDDO
           ELSE IF (IDK_YT(K).EQ.1) THEN
              DO I=2,N
                 NN=NN+1
                 YY(NN)=X(I,IDX_YT(K))
            ENDDO
           ENDIF              
        ENDDO
        NN_YT=NN
        ELSE   ! (1)
      NN=0
      DO K=1,12
         IF (IDK_YT(K).EQ.0) THEN
              DO I=1,N
                 NN=NN+1
                 YY(NN)=X(I,IDX_YT(K))
              ENDDO
           ELSE IF (IDK_YT(K).EQ.1) THEN
              DO I=1,N
                 NN=NN+1
                 YY(NN)=X(I,IDX_YT(K))
            ENDDO
           ENDIF              
        ENDDO
        NN_YT=NN
        END IF ! (1)
C
C...READ PREDICTOR-DATA 
C
      DO 50 L=1,MVAR     

        OPEN (41,FILE=FNAME_NT(L),STATUS='OLD')
        M=MMAX
        READ(41,*)
        READ(41,*) !N, YEAR1, YEAR2
        READ(41,*)
        DO I=1,N
           READ(41,*) (X(I,J),J=0,M)
        ENDDO

        IF (SHIFT.EQ.1) THEN ! (2)
      NN=0
      DO K=1,12
         IF (IDK_NT(K).EQ.0) THEN
              DO I=1,N-1
                 NN=NN+1
                 XX(NN,L)=X(I,IDX_NT(K))
              ENDDO
           ELSE IF (IDK_NT(K).EQ.1) THEN
              DO I=2,N
                 NN=NN+1
                 XX(NN,L)=X(I,IDX_NT(K))
            ENDDO
           ENDIF              
        ENDDO
        NN_NT=NN
        ELSE   ! (2)
      NN=0
      DO K=1,12
         IF (IDK_NT(K).EQ.0) THEN
              DO I=1,N
                 NN=NN+1
                 XX(NN,L)=X(I,IDX_NT(K))
              ENDDO
           ELSE IF (IDK_NT(K).EQ.1) THEN
              DO I=1,N
                 NN=NN+1
                 XX(NN,L)=X(I,IDX_NT(K))
            ENDDO
           ENDIF              
        ENDDO
        NN_NT=NN
        END IF ! (2)
        CLOSE(41)
  50  CONTINUE
C
C   TACH NHOM
C
      N1=0
        N2=0
        N3=0
        DO I=1,NN_NT
           IF (YY(I).EQ.1.0) THEN
              N1=N1+1
              DO L=1,MVAR
                 XG(N1,L,1) = XX(I,L)
                  ENDDO
           ELSE IF (YY(I).EQ.2.0) THEN
              N2=N2+1
              DO L=1,MVAR
             XG(N2,L,2) = XX(I,L)
              ENDDO
           ELSE IF (YY(I).EQ.3.0) THEN
              N3=N3+1
              DO L=1,MVAR
                 XG(N3,L,3) = XX(I,L)
              ENDDO
         ENDIF
        ENDDO
C
        OPEN (43,FILE='TESTOUT.TXT', STATUS='UNKNOWN')
      WRITE(43,*) KMAX,MVAR,N1,N2,N3
C... LOP 1
           DO I=1,N1
              WRITE(43,'(8F7.2)') (XG(I,L,1),L=1,MVAR)
           ENDDO

C... LOP 2
           DO I=1,N2
              WRITE(43,'(8F7.2)') (XG(I,L,2),L=1,MVAR)
           ENDDO

C... LOP 3
           DO I=1,N3
              WRITE(43,'(8F7.2)') (XG(I,L,3),L=1,MVAR)
           ENDDO
        CLOSE (43)
        RETURN
        END ! TACH SO LIEU

!###########################################
        SUBROUTINE MAKECODE(FNAME)
        PARAMETER (NMAX=120,MMAX=12)
        CHARACTER FNAME*20, ST1*100, ST2*100
        REAL XX(NMAX,0:MMAX), X(NMAX)
        REAL Q33,Q66
        INTEGER N,I,J,YEAR1, YEAR2

        OPEN (15,FILE=FNAME,STATUS='OLD')

        READ (15,'(A100)') ST1
        READ (15,*) N,YEAR1,YEAR2
        READ (15,'(A100)') ST2
        DO I=1,N
           READ (15,*) (XX(I,J),J=0,MMAX)
        ENDDO
        CLOSE(15)

        DO J=1,MMAX
           DO I=1,N
              X(I)=XX(I,J)
           ENDDO
           CALL QUANTILE (X,N,Q33,Q66)
           DO I=1,N
              IF (XX(I,J).LE.Q33) THEN
                 XX(I,J)=1.0
              ELSEIF (XX(I,J).LE.Q66) THEN
                 XX(I,J)=2.0
              ELSE
                 XX(I,J)=3.0
              ENDIF
           ENDDO
        ENDDO

        OPEN (15,FILE='DRAFT.TXT',STATUS='UNKNOWN')

        WRITE (15,'(A100)') ST1
        WRITE (15,'(3I6)') N, YEAR1,YEAR2
        WRITE (15,'(A100)') ST2
        DO I=1,N
           WRITE (15,'(13F7.0)') (XX(I,J),J=0,MMAX)
        ENDDO
        CLOSE(15)

        RETURN
        END

!###########################################
      SUBROUTINE QUANTILE (X,N,Q33,Q66)
C
C  TINH PHAN VI 33% VA 66% CUA CHUOI X
C
      REAL X(N), W(N), Q33,Q66
        W = X
      DO I=1,N-1
        DO J=I+1,N
           IF (W(I).GT.W(J)) THEN
              TMP=W(I)
              W(I)=W(J)
              W(J)=TMP
           ENDIF
        ENDDO
        ENDDO
        M = N/3
      K = MOD (N,M)
        IF (K.EQ.0) THEN
           Q33=(W(M)+W(M+1))/2.
         Q66=(W(2*M)+W(2*M+1))/2.
        ELSE IF (K.EQ.1) THEN
           Q33=(W(M)+W(M+1))/2.
         Q66=(W(2*M+1)+W(2*M+2))/2.
        ELSE IF (K.EQ.2) THEN
           Q33= W(M+1)
         Q66= W(2*M+2)
      ENDIF
        RETURN
        END

!###########################################

      SUBROUTINE CHON_TINH
      character *12 fsl,fra,FTEMP
        INTEGER NGROUP, NVAR
        PARAMETER (NGROUP=3, NVAR=8, NMAX=500) ! SO NHOM, SO BIEN, 
                                               ! SO MAU CUC DAI MOI NHOM
        PARAMETER (MVAR1=8,NTOHOP=400, NNHOM=3)
        REAL PTPL(NNHOM,MVAR1+1), DCX
      REAL PTR(NTOHOP,NNHOM,MVAR1+1), ACCURATE(NTOHOP) ! LUU CAC PHUONG

        INTEGER NT(NGROUP)
        REAL X(NGROUP*NMAX, NVAR),Y(NGROUP*NMAX, NVAR)

        COMMON /STORE_KQ/ PTR, ACCURATE, IPT
        COMMON /ACCU/ PTPL,DCX
        COMMON /FNAME/ fsl,fra,FTEMP
        INTEGER I1,I2,I3,I4,I5,I6,I7,I8
        COMMON /VARIDX/ I1,I2,I3,I4,I5,I6,I7,I8

      IPT=0   ! DEM SO PHUONG TRINH
        I1=0
        I2=0
        I3=0
        I4=0
        I5=0
        I6=0
        I7=0
        I8=0

        OPEN (31, FILE='TESTOUT.TXT',STATUS='OLD')
        READ (31,*) TMP, TMP, (NT(I),I=1,NGROUP)
        N=0
        DO I=1,NGROUP
           N=N+NT(I)  ! TONG DUNG LUONG MAU CAC NHOM
        ENDDO

        DO I=1,N
           READ (31,*) (X(I,J),J=1,NVAR) ! SO LIEU CAC BIEN
        ENDDO
        CLOSE(31)

! HAI BIEN
        PRINT*,' CALCULATING 2 VARIABLES....'
      MVARS=2
      DO I1=1,NVAR-1
           DO I2=I1+1,NVAR
        OPEN (32, FILE='TEMP.TXT',STATUS='UNKNOWN')
            DO I=1,N
                 Y(I,1)=X(I,I1)
                 Y(I,2)=X(I,I2)
            ENDDO
              WRITE (32,'(20I5)') NGROUP,MVARS,(NT(J),J=1,NGROUP)
              DO I=1,N
                 WRITE (32,'(20F10.2)') (Y(I,J),J=1,MVARS)
              ENDDO
        CLOSE(32)
              CALL PLMAIN (MVARS,FTEMP)
              IPT=IPT+1
              CALL STORE (IPT,MVARS)
           ENDDO
        ENDDO
        
! BA BIEN
        PRINT*,' CALCULATING 3 VARIABLES....'
      MVARS=3
      DO I1=1,NVAR-2
           DO I2=I1+1,NVAR-1
              DO I3=I2+1,NVAR
        OPEN (32, FILE='TEMP.TXT',STATUS='UNKNOWN')
            DO I=1,N
                 Y(I,1)=X(I,I1)
                 Y(I,2)=X(I,I2)
                 Y(I,3)=X(I,I3)
            ENDDO
              WRITE (32,'(20I5)') NGROUP,MVARS,(NT(J),J=1,NGROUP)
              DO I=1,N
                 WRITE (32,'(20F10.2)') (Y(I,J),J=1,MVARS)
              ENDDO
      CLOSE(32)
              CALL PLMAIN (MVARS,FTEMP)
              IPT=IPT+1
              CALL STORE (IPT,MVARS)
              ENDDO
           ENDDO
        ENDDO

! BON BIEN
        PRINT*,' CALCULATING 4 VARIABLES....'
      MVARS=4
      DO I1=1,NVAR-3
           DO I2=I1+1,NVAR-2
              DO I3=I2+1,NVAR-1
                 DO I4=I3+1,NVAR
        OPEN (32, FILE='TEMP.TXT',STATUS='UNKNOWN')
            DO I=1,N
                 Y(I,1)=X(I,I1)
                 Y(I,2)=X(I,I2)
                 Y(I,3)=X(I,I3)
                 Y(I,4)=X(I,I4)
            ENDDO
              WRITE (32,'(20I5)') NGROUP,MVARS,(NT(J),J=1,NGROUP)
              DO I=1,N
                 WRITE (32,'(20F10.2)') (Y(I,J),J=1,MVARS)
              ENDDO
        CLOSE(32)
              CALL PLMAIN (MVARS,FTEMP)
              IPT=IPT+1
              CALL STORE (IPT,MVARS)
                 ENDDO
              ENDDO
           ENDDO
        ENDDO

! NAM BIEN
        PRINT*,' CALCULATING 5 VARIABLES....'
      MVARS=5
      DO I1=1,NVAR-4
           DO I2=I1+1,NVAR-3
              DO I3=I2+1,NVAR-2
                 DO I4=I3+1,NVAR-1
                    DO I5=I4+1,NVAR
        OPEN (32, FILE='TEMP.TXT',STATUS='UNKNOWN')
            DO I=1,N
                 Y(I,1)=X(I,I1)
                 Y(I,2)=X(I,I2)
                 Y(I,3)=X(I,I3)
                 Y(I,4)=X(I,I4)
                 Y(I,5)=X(I,I5)
            ENDDO
              WRITE (32,'(20I5)') NGROUP,MVARS,(NT(J),J=1,NGROUP)
              DO I=1,N
                 WRITE (32,'(20F10.2)') (Y(I,J),J=1,MVARS)
              ENDDO
        CLOSE(32)
              CALL PLMAIN (MVARS,FTEMP)
              IPT=IPT+1
              CALL STORE (IPT,MVARS)
                    ENDDO
                 ENDDO
              ENDDO
           ENDDO
        ENDDO

! SAU BIEN
        PRINT*,' CALCULATING 6 VARIABLES....'
      MVARS=6
      DO I1=1,NVAR-5
           DO I2=I1+1,NVAR-4
              DO I3=I2+1,NVAR-3
                 DO I4=I3+1,NVAR-2
                    DO I5=I4+1,NVAR-1
                       DO I6=I5+1,NVAR
        OPEN (32, FILE='TEMP.TXT',STATUS='UNKNOWN')
            DO I=1,N
                 Y(I,1)=X(I,I1)
                 Y(I,2)=X(I,I2)
                 Y(I,3)=X(I,I3)
                 Y(I,4)=X(I,I4)
                 Y(I,5)=X(I,I5)
                 Y(I,6)=X(I,I6)
            ENDDO
              WRITE (32,'(20I5)') NGROUP,MVARS,(NT(J),J=1,NGROUP)
              DO I=1,N
                 WRITE (32,'(20F10.2)') (Y(I,J),J=1,MVARS)
              ENDDO
        CLOSE(32)
              CALL PLMAIN (MVARS,FTEMP)
              IPT=IPT+1
              CALL STORE (IPT,MVARS)
                       ENDDO
                    ENDDO
                 ENDDO
              ENDDO
           ENDDO
        ENDDO

! BAY BIEN
        PRINT*,' CALCULATING 7 VARIABLES....'
      MVARS=7
      DO I1=1,NVAR-6
           DO I2=I1+1,NVAR-5
               DO I3=I2+1,NVAR-4
                 DO I4=I3+1,NVAR-3
                    DO I5=I4+1,NVAR-2
                       DO I6=I5+1,NVAR-1
                          DO I7=I6+1,NVAR
        OPEN (32, FILE='TEMP.TXT',STATUS='UNKNOWN')
            DO I=1,N
                 Y(I,1)=X(I,I1)
                 Y(I,2)=X(I,I2)
                 Y(I,3)=X(I,I3)
                 Y(I,4)=X(I,I4)
                 Y(I,5)=X(I,I5)
                 Y(I,6)=X(I,I6)
                 Y(I,7)=X(I,I7)
            ENDDO
              WRITE (32,'(20I5)') NGROUP,MVARS,(NT(J),J=1,NGROUP)
              DO I=1,N
                 WRITE (32,'(20F10.2)') (Y(I,J),J=1,MVARS)
              ENDDO
        CLOSE(32)
              CALL PLMAIN (MVARS,FTEMP)
              IPT=IPT+1
              CALL STORE (IPT,MVARS)
                          ENDDO
                       ENDDO
                    ENDDO
                 ENDDO
              ENDDO
           ENDDO
        ENDDO
        RETURN
        END
!##################################################

      SUBROUTINE STORE (IPT,MVAR)
        PARAMETER (MVAR1=8,NTOHOP=400, NNHOM=3)
        REAL PTPL(NNHOM,MVAR1+1), DCX
      REAL PTR(NTOHOP,NNHOM,MVAR1+1), ACCURATE(NTOHOP) ! LUU CAC PHUONG
        INTEGER I1,I2,I3,I4,I5,I6,I7,I8
        INTEGER IDXVAR(MVAR1+1)

        COMMON /VARIDX/ I1,I2,I3,I4,I5,I6,I7,I8
        COMMON /STORE_KQ/ PTR, ACCURATE
        COMMON /ACCU/ PTPL,DCX

        IDXVAR(1)= 1
        IDXVAR(2)=I1 + 1
        IDXVAR(3)=I2 + 1
        IDXVAR(4)=I3 + 1
        IDXVAR(5)=I4 + 1
        IDXVAR(6)=I5 + 1
        IDXVAR(7)=I6 + 1
        IDXVAR(8)=I7 + 1
        IDXVAR(9)=I8 + 1
!       PRINT '(10I5)',(IDXVAR(I),I=1,9)
!       PAUSE

        DO I=1,NNHOM
!          DO J=1,MVAR+1
           DO J=1,MVAR1+1
              IF (J.EQ.1) THEN 
                     PTR(IPT,I,IDXVAR(J))=PTPL(I,J)
              ELSE
                 IF (IDXVAR(J).NE.1) PTR(IPT,I,IDXVAR(J))=PTPL(I,J)  
                                  ! LUU CAC PT TRA VE CHUONG TRINH CHINH
              ENDIF
           ENDDO
        ENDDO
        ACCURATE(IPT)=DCX
      RETURN
        END
C
!###################################################
      SUBROUTINE PLMAIN (MVAR,FTEMP)
                 ! CHUONG TRINH TINH HAM PHAN LOP
        PARAMETER (MVAR1=8,NTOHOP=400, NNHOM=3)
        PARAMETER (NLOP=10)
      dimension N(10),CMEAN(10),XBAR(100),C(110),D(100)
      DIMENSION P(500),X(5000)
        integer LG(500)
      character *12 fsl,FTEMP
        DIMENSION ND(NLOP),ND1(NLOP),ND2(NLOP),ND3(NLOP)
        REAL PP(NLOP,NLOP), P1(NLOP),P2(NLOP)   ! LUU KET QUA DANG GIA
        REAL PTPL(NNHOM,MVAR1+1), DCX

        COMMON /ACCU/ PTPL,DCX
      COMMON /FNAME/ fsl
C
C   Kich thuoc duoc mo ta nhu sau
C   Cua N : >= so nhom K
C   Cua CMEAN : >= so bien M
C   Cua XBAR  : >= tich so M x K
C   Cua C : >= tich (M+1) x K
C   Cua D : >= tich M x M
C   cua P,LG : >= kich thuoc chung cua tap mau K nhom, T
C                (T=N(1)+N(2)+...+N(k))
C   Cua X : >= tich T x M
C       

      open(1,File=FTEMP,status='OLD')

  100 read(1,*) K,M,(N(i),i=1,K)
      write(3,2) K,M
    2 format('Chuong trinh Ham phan lop'/'So nhom',I5/'So bien',I5/)

      L=0
      do 130 i=1,k
         n1=n(i)
         do 120 j=1,n1
            read(1,*)(cmean(ij),ij=1,M)
            L=L+1
            n2=L-n1
            do 120 ij=1,M
               n2=n2+n1
  120    x(n2)=cmean(ij)
  130 L=n2

      call DMATX(k,m,n,x,xbar,d,cmean)
      L=0
      do 150 i=1,k
         do 140 j=1,m
            L=L+1
  140    cmean(j)=xbar(L)
  150 continue

      do 170 i=1,m
         L=i-m
         do 160 j=1,m
            L=L+m
  160    cmean(j)=D(L)

  170 continue

      Call MINV(D,M,DET,CMEAN,C)
      Call DISCR(K,M,N,X,XBAR,D,CMEAN,V,C,P,LG)

    8 format(' Hang ',i3/(8f15.5))
    9 format(' Trung binh chung '/(8f15.5))
   10 format(' Hang so D-Mahalanov ',f15.5)

        WRITE (3,*)' TRUONG HOP ', MVAR, ' BIEN'
!###############################################
      n1=1
      n2=M+1
      do 180 i=1,k
         write(3,11) I,(C(J),J=N1,N2)
           I1=0
           DO J=N1,N2
              I1=I1+1
              PTPL(I,I1)=C(J)   ! LUU PHUONG TRINH PHAN LOP
           ENDDO
   11    FORMAT(' Ham phan lop ',i3/6x,'Hang so    Cac he so'/9f12.3)
         n1=n1+(M+1)
  180 n2=n2+(M+1)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ND=0
        ND1=0
        ND2=0
        ND3=0
        PP = 0.0 ! GAN BANG 0 CHO MANG KET QUA DANH GIA
        P1 = 0.0
        P2 = 0.0
        TONG = 0.0
         
      n1=1
      n2=N(1)
        do 210 i=1,k
                
!##############################################
                DO II=1,K
             DO J=N1,N2
                  IF (LG(J).EQ.II) then 
                   PP(I,II)=PP(I,II)+1
                  ENDIF
               ENDDO
            ENDDO
! ######################################
                L=0
                do 190 j=n1,n2
                   L=L+1
                   IF (LG(J).EQ.1) then 
                      ND1(I)=ND1(I)+1
               else if (LG(J).EQ.2) then 
                      ND2(I)=ND2(I)+1 
               else 
                      ND3(I)=ND3(I)+1
               endif
  190       CONTINUE !!         write(3,14) L,LG(j)

   14           format(1x,i4,2x,i6)

                if (i-k) 200,300,300
  200           n1=n1+N(i)
                n2=n2+N(i+1)
  210   continue
  300 continue
!###########################################
          DCX=0.0
            DO J=1,K    ! CONG HANG DUOI CUNG (CONG THEO COT)
               DO II=1,K
                  P1(J)=P1(J)+PP(II,J)
                  TONG=TONG+PP(II,J)
               ENDDO
            ENDDO
            DO II=1,K    ! CONG COT BEN PHAI NHAT (CONG THEO HANG)
               DO J=1,K
                  P2(II)=P2(II)+PP(II,J)
               ENDDO
            ENDDO
            DO II=1,K
               P1(II)=P1(II)/TONG*100
               P2(II)=P2(II)/TONG*100
               DO J=1,K
                  PP(II,J)=PP(II,J)/TONG*100
                  IF (II.EQ.J) DCX=DCX+PP(II,J)
               ENDDO
            ENDDO
!##############################################
      DO I=1,K
          WRITE (3,'(10F10.1)') (PP(I,II),II=1,K), P2(I)
          ENDDO
          WRITE (3,'(10F10.1)') (P1(I),I=1,K), TONG
          WRITE (3,'("    DO CHINH XAC CHUNG = ",F10.2)') DCX
          WRITE (3,*)'******************************'
          CLOSE(1)
      end


C      Chuong trinh con trong bo chuong trinh
C      phan tich phan biet (phan lop)
C
      SUBROUTINE DMATX(K,M,N,X,XBAR,D,CMEAN)
      dimension N(1),X(1),XBAR(1),D(1),CMEAN(1)
      MM=M*M
      do 100 i=1,MM
  100 D(i)=0
      N4=0
      L=0
      LM=0
      do 160 NG=1,K
      N1=N(NG)
      FN=N1
      do 130 j=1,M
      LM=LM+1
      XBAR(LM)=0
      do 120 i=1,N1

      L=L+1
  120 XBAR(LM)=XBAR(LM)+X(L)
  130 XBAR(LM)=XBAR(LM)/FN
      LMEAN=LM-M
      do 150 i=1,N1
      LL=N4+i-N1
      do 140 j=1,M
      LL=LL+N1
      N2=LMEAN+j
  140 CMEAN(j)=X(LL)-XBAR(N2)
      LL=0
      do 150 j=1,M
      do 150 jj=1,M
      LL=LL+1
  150 D(LL)=D(LL)+CMEAN(J)*CMEAN(JJ)
  160 N4=N4+N1*M
      LL=-K
      DO 170 I=1,K
  170 LL=LL+N(I)
      FN=LL
      do 180 I=1,MM
  180 D(I)=D(I)/FN
      RETURN
      END


      SUBROUTINE MINV(A,N,D,L,M)
      DIMENSION A(1),L(1),M(1)
      D=1.
      NK=-N
      DO 80 K=1,N
      NK=NK+N
      L(K)=K
      M(K)=K
      KK=NK+K
      BIGA=A(KK)
      DO 20 J=K,N
      IZ=N*(J-1)
      DO 20 I=K,N
      IJ=IZ+I
   10 IF(ABS(BIGA)-ABS(A(IJ))) 15,20,20
   15 BIGA=A(IJ)
      L(K)=I
      M(K)=J
   20 CONTINUE
      J=L(K)
      IF(J-K)35,35,25
   25 KI=K-N
      DO 30 I=1,N
      KI=KI+N
      HOLD=-A(KI)
      JI=KI-K+J
      A(KI)=A(JI)
   30 A(JI)=HOLD
   35 I=M(K)
      IF(I-K)45,45,38
   38 JP=N*(I-1)
      DO 40 J=1,N
      JK=NK+J
      JI=JP+J
      HOLD=-A(JK)
      A(JK)=A(JI)
   40 A(JI)=HOLD
   45 IF(BIGA)48,46,48
   46 D=0.
      RETURN
   48 DO 55 I=1,N
      IF(I-K)50,55,50
   50 IK=NK+I
      A(IK)=A(IK)/(-BIGA)
   55 CONTINUE
      DO 65 I=1,N
      IK=NK+I
      HOLD=A(IK)
      IJ=I-N
      DO 65 J=1,N
      IJ=IJ+N
      IF(I-K)60,65,60
   60 IF(J-K)62,65,62
   62 KJ=IJ-I+K
      A(IJ)=HOLD*A(KJ)+A(IJ)
   65 CONTINUE
      KJ=K-N
      DO 75 J=1,N
      KJ=KJ+N
      IF(J-K)70,75,70
   70 A(KJ)=A(KJ)/BIGA
   75 CONTINUE
      D=D*BIGA
      A(KK)=1./BIGA
   80 CONTINUE
      K=N
  100 K=(K-1)
      IF(K)150,150,105
  105 I=L(K)
      IF(I-K)120,120,108
  108 JQ=N*(K-1)
      JR=N*(I-1)
      DO 110 J=1,N
      JK=JQ+J
      HOLD=A(JK)
      JI=JR+J
      A(JK)=-A(JI)
  110 A(JI)=HOLD
  120 J=M(K)
      IF(J-K)100,100,125
  125 KI=K-N
      DO 130 I=1,N
      KI=KI+N
      HOLD=A(KI)
      JI=KI-K+J
      A(KI)=-A(JI)
  130 A(JI)=HOLD
      GOTO 100
  150 RETURN
        END

C     Chuong trinh trong bo chuong trinh phan lop

      SUBROUTINE DISCR(K,M,N,X,XBAR,D,CMEAN,V,C,P,LG)
      dimension N(1),X(1),XBAR(1),D(1),CMEAN(1),C(1),P(1),LG(1)
      N1=N(1)
      do 100 i=2,k
  100 N1=N1+N(I)
      FN1=N1
      do 110 i=1,k
  110 P(i)=N(I)
      do 130 i=1,M
      CMEAN(I)=0
      N1=I-M
      do 120 j=1,k
      N1=N1+M
  120 CMEAN(I)=CMEAN(I)+P(J)*XBAR(N1)
  130 CMEAN(I)=CMEAN(I)/FN1
      L=0
      do 140 i=1,k
      do 140 j=1,M
      L=L+1
  140 C(L)=XBAR(L)-CMEAN(J)
      V=0
      L=0
      do 160 j=1,M
      do 160 i=1,M
      N1=I-M
      N2=J-M
      Sum=0
      do 150 ij=1,k
      N1=N1+M
      N2=N2+M
  150 Sum=Sum+P(IJ)*C(N1)*C(N2)
      L=L+1
  160 V=V+D(L)*Sum
      N2=0
      do 190 KA=1,k
      do 170 i=1,M
      N2=N2+1
  170 P(I)=XBAR(N2)
      IQ=(M+1)*(KA-1)+1
      Sum=0
      do 180 j=1,M
      N1=J-M
      do 180 L=1,M
      N1=N1+M
  180 Sum=Sum+D(N1)*P(J)*P(L)
      C(IQ)=-(Sum/2)
      do 190 i=1,M
      N1=I-M
      IQ=IQ+1
      C(IQ)=0
      do 190 j=1,M
      N1=N1+M
  190 C(IQ)=C(IQ)+D(N1)*P(J)
      LBASE=0
      N1=0
      do 270 KG=1,K
      NN=N(KG)
      do 260 i=1,NN
      L=I-NN+LBASE
      do 200 j=1,M
      L=L+NN
  200 D(J)=X(L)
      N2=0
      do 220 KA=1,K
      N2=N2+1
      Sum=C(N2)
      do 210 j=1,M
      N2=N2+1
  210 Sum=Sum+C(N2)*D(J)
  220 XBAR(KA)=Sum
      L=1
      Sum=XBAR(1)
      do 240 j=2,k
      IF(Sum-XBAR(j)) 230,240,240
  230 L=J
      Sum=XBAR(J)
  240 CONTINUE
      PL=0
      do 250 j=1,k
  250 PL=PL+EXP(XBAR(J)-Sum)
      N1=N1+1
      LG(N1)=L
  260 P(N1)=1/PL
  270 LBASE=LBASE+NN*M
      RETURN
      END
!######################################################