$DEBUG
      PROGRAM MAINPRG
      INCLUDE 'PLOP.INC'
!       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
        INTEGER NTERM, NMONTH
        DIMENSION NDLG(NNHOM)                   ! DUNG LUONG MAU CAC
NHOM
        
      REAL PTR(NTOHOP,NNHOM,MVAR1+1), ACCURATE(NTOHOP) ! LUU CAC PHUONG
TRINH
                                                                                                  ! VA
                                                                                                  ! DO
                                                                                                  ! CHINH
                                                                                                  ! XAC
        CHARACTER HSPTR(NTOHOP,NNHOM,MVAR1+1)*10, VARNAME(MVARMAX)*10
        
      INTEGER IPT, N_YEUTO, N_VUNG
        CHARACTER *70 FNAME_NT(MVARMAX), FNAME_YT, FN_YT(NVUNG)
        CHARACTER FOUT*70,ST*30,ST1*70, FOUT1*70, F_DCX*70
        CHARACTER *70 F_OUT(NVUNG), F_OUT1(NVUNG)
        CHARACTER *70 DIR_YEUTO(NYEUTO),DIR_OUT(NYEUTO)

      character *70 fsl,fra,FTEMP
        CHARACTER ST2*100,TIEUDE*100
        LOGICAL QUANTILE
        REAL PP0(NTOHOP,NNHOM,NNHOM),P10(NTOHOP,NNHOM),P20(NTOHOP,NNHOM)
        INTEGER VUNG, YEUTO

        NAMELIST  /FN_IN_OUT/ FSL, FTEMP, FRA 
        1          /CONTROL/   NMONTH, QUANTILE,        N_VUNG, N_YEUTO,
        1                      DIR_YEUTO,DIR_OUT
        2          /FNAMENT/   FNAME_NT
        3          /FNAMEYT/   FN_YT, F_OUT, FOUT1, F_DCX, TIEUDE
        4          /TENBIEN/   VARNAME

        COMMON /TABLE/ PP0,P10,P20
        COMMON /STORE_KQ/ PTR, ACCURATE, IPT
        COMMON /FNAME/   fsl,fra,FTEMP
        COMMON /CTRL/    IDX_YT,IDK_YT,IDX_NT,IDK_NT,NTERM, NMONTH
        COMMON /FNAM_NT/ FNAME_NT
        COMMON /FNAM_YT/ FNAME_YT, FOUT, TIEUDE
        COMMON /DLMAU/ NDLG

      OPEN (21,FILE='INPUT_DISR2.TXT',STATUS='OLD')
        READ (21,FN_IN_OUT)
!       WRITE (*,FN_IN_OUT)
        READ (21,CONTROL)
!       WRITE (*,CONTROL)
      READ (21,FNAMENT  )
!      WRITE (*,FNAMENT  )
      READ (21,FNAMEYT  )
!      WRITE (*,FNAMEYT  )
        READ (21,TENBIEN  )
!       WRITE (*,TENBIEN  )

        CLOSE (21)
        
        OPEN (56, FILE=F_DCX, STATUS='UNKNOWN')
      OPEN (55,FILE=FOUT1,STATUS='UNKNOWN',ACCESS='APPEND')

        DO 1001 YEUTO = 1, N_YEUTO
           LENDIRYT = INDEX(DIR_YEUTO(YEUTO),' ')-1
           LENDIROUT= INDEX(DIR_OUT(YEUTO),' ')-1
        
      DO 1001 VUNG  = 1, N_VUNG

           FNAME_YT = DIR_YEUTO(YEUTO)(1:LENDIRYT)//FN_YT(VUNG)
           FOUT     = DIR_OUT(YEUTO)(1:LENDIROUT)//F_OUT(VUNG)
!       PRINT*,'VUNG TINH : ',FNAME_YT
!       PRINT*,'YEU TO : ',YEUTO
!       PRINT*,'OUT FILE: ',FOUT
!       PRINT*,'TOTAL OUT:',FOUT1
!       PAUSE
        
        PTR=-999.0
        HSPTR='      ----'

        ST1=FOUT
        ST2=TIEUDE
        LENTD=INDEX(TIEUDE,'#')-1
        PRINT*,LENTD
      open (3,File=FRA,status='unknown')
        OPEN (54,FILE=FOUT,STATUS='UNKNOWN')

C..LOOP
        DO 1000 ICYCLE=1, 12
        DO 1000 NTERM =1,3

        PRINT*,'=========== ', ICYCLE,' =========='
        WRITE(ST,'("_",I2.2,"_",I1.1,"_",I1.1)') ICYCLE,NMONTH,NTERM

      CALL MAKEIDX (ICYCLE)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++
        CALL TACHSL (QUANTILE)    !...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

        LENST1=INDEX(ST1,' ')-1
        LENST =INDEX(ST ,' ')-1
        FOUT=ST1(1:LENST1)//ST(1:LENST)//'.TXT'
        
!!      OPEN (55,FILE=FOUT,STATUS='UNKNOWN')
        
      ST2=TIEUDE(1:LENTD)//' ('//ST(1:LENST)//')'
!!      WRITE (55,'(A100)') ST2
      WRITE (54,'(A100)') ST2
!!      WRITE (55,'("     ","   HESO_Co",10A10)') (VARNAME(I),I=1,MVAR1)
      DO I=1,IPT
!C         WRITE(54,'("PTRINH ",I3)') I
!C         WRITE(54,'("DCXAC: ",F10.1)') ACCURATE(I)
!C         DO J=1,NNHOM
!C            WRITE (54,'("HAM",I2,10A10)') J,(HSPTR(I,J,K),K=1,MVAR1+1)
!C         ENDDO

           IF (ACCURATE(I).GT.DCXMAX) THEN
              DCXMAX=ACCURATE(I)
              IPTMAX=I
           ENDIF
        ENDDO

!!      DO J=1,NNHOM
!!         WRITE (55,'("HAM",I2,10A10)')
!J,(HSPTR(IPTMAX,J,K),K=1,MVAR1+1)
!!      ENDDO

        DO K=1, MVAR1 + 1
           IF (K.EQ.1) THEN
              WRITE (55,'(4I4,"  HE_S0_A0",3A10)')
VUNG,YEUTO,ICYCLE,NTERM
        &            , (HSPTR(IPTMAX,J,K),J=1, NNHOM)
           ELSE
              WRITE (55,'(4I4,4A10)')
VUNG,YEUTO,ICYCLE,NTERM,VARNAME(K-1)
        &            , (HSPTR(IPTMAX,J,K),J=1, NNHOM)
           ENDIF
        ENDDO

!!
        WRITE (54,*) '======================================'
        WRITE (54,'(" P.TRINH THU ",I3," - DO C/XAC MAX = ",F10.4)') 
        &              IPTMAX, DCXMAX
        WRITE (54,*)
        WRITE (54,'("     ",10I10)') (NDLG(IDL),IDL=1,NNHOM)
        WRITE (54,*)
      DO J=1,NNHOM
      WRITE (54,'("   ",I2,10F10.1)') J,(PP0(IPTMAX,J,K),K=1,NNHOM),
        &                                      P20(IPTMAX,J)
        ENDDO

        WRITE (54,*)
        WRITE (54,'("     ",10F10.1)') (P10(IPTMAX,J),J=1,NNHOM)

C...WRITE ACCURACY TO FILE DCX.TXT
C*********************************
        WRITE (56, '(4I4,F10.4)') VUNG,YEUTO,ICYCLE,NTERM, DCXMAX

!!
        PRINT*, 'TRUONG HOP ',ICYCLE, ' - DLG MAU=', 
        &                          (NDLG(IDL),IDL=1,NNHOM)
        PRINT '(" P.TRINH THU ",I3," - DO C/XAC MAX = ",F10.4)', 
        &              IPTMAX, DCXMAX

 1000 CONTINUE
 1001 CONTINUE

      END ! MAIN PROGRAM
!####################################################
        SUBROUTINE MAKEIDX (I)
        INTEGER IDX_YT(12), IDK_YT(12), IDX_NT(12), IDK_NT(12)
        INTEGER NTERM, NMONTH
        COMMON /CTRL/    IDX_YT,IDK_YT,IDX_NT,IDK_NT,NTERM, NMONTH
           IDX_YT=0
           IDX_NT=0
           IDK_YT=0
           IDK_NT=0
           DO J=1,NMONTH
              IF (J.EQ.1) THEN 
                     IDX_YT(J) = I
              ELSE
                 IDX_YT(J) = IDX_YT(J-1) + 1
                 IF (IDX_YT(J).GT.12) IDX_YT(J)=1
              ENDIF

              IF (J.EQ.1) THEN
                 IDX_NT(J) = IDX_YT(J)-NTERM-1
                 IF (IDX_NT(J).LE.0) IDX_NT(J)=IDX_NT(J)+12
                 IF (IDX_NT(J).GT.12) IDX_NT(J)=1
              ELSE
                 IDX_NT(J) = IDX_NT(J-1) + 1
                 IF (IDX_NT(J).GT.12) IDX_NT(J)=1
              ENDIF

              IF (IDX_YT(J).LT.IDX_NT(J)) THEN
                 IDK_YT(J)=1
                 IDK_NT(J)=1
              ENDIF
           ENDDO
        RETURN
        END

!####################################################
      SUBROUTINE TACHSL (QUANTILE)
!       PARAMETER (MVAR1=8,NTOHOP=400, NNHOM=3, NMAX=200, MMAX=12)
      INCLUDE 'PLOP.INC'

        PARAMETER (KMAX=NNHOM,MVAR=MVAR1)
        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 NN,YEAR1,YEAR2
        CHARACTER *70 FNAME_NT(MVAR), FNAME_YT
        LOGICAL QUANTILE

        COMMON /CTRL/    IDX_YT,IDK_YT,IDX_NT,IDK_NT,NTERM, NMONTH
        COMMON /FNAM_NT/ FNAME_NT
        COMMON /FNAM_YT/ FNAME_YT

C...MAKE CODE
        CALL MAKECODE(FNAME_YT,QUANTILE)

C...READ PREDICTANT-DATA 

        N = 41 ! DUNG LUONG MAU CHUNG
C
        OPEN (41,FILE='DRAFT.TXT',STATUS='OLD')

        READ(41,*)
        READ(41,*) ! N, YEAR1, YEAR2
        READ(41,*)

        M=MMAX
        DO I=1,N
           READ(41,*) (X(I,J),J=0,M)
        ENDDO
        CLOSE(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=2,N
                 NN=NN+1
                 YY(NN)=X(I,IDX_YT(K))
            ENDDO
           ENDIF              
        ENDDO
        NN_YT=NN
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

      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-1
                 NN=NN+1
                 XX(NN,L)=X(I,IDX_NT(K))
            ENDDO
           ENDIF              
        ENDDO
        NN_NT=NN
        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,QUANTILE)

      INCLUDE 'PLOP.INC'
!       PARAMETER (MVAR1=8,NTOHOP=400, NNHOM=3, NMAX=200, MMAX=12)
!       PARAMETER (NMAX=120,MMAX=12)
        CHARACTER FNAME*70, ST1*100, ST2*100
        REAL XX(NMAX,0:MMAX), X(NMAX)
        REAL Q33,Q66
        INTEGER N,I,J,YEAR1, YEAR2
        LOGICAL QUANTILE

      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
           X=-999.0
           DO I=1,N
              X(I)=XX(I,J)
           ENDDO
           IF (QUANTILE) THEN
              CALL QUANTILES (X,N,Q33,Q66)
           ELSE
              CALL STDEV     (X,N,Q33,Q66)
           ENDIF
           DO I=1,N
              IF (XX(I,J).LE.Q33) THEN
                 XX(I,J)=1.0
              ELSEIF (XX(I,J).LT.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 QUANTILES (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 STDEV (X,N,Q33,Q66)
C
C  TINH NGUONG CHIA THEO STANDARD DIVIATION
C
      REAL X(N), W(N), Q33,Q66

        W = X
        SUM1=0.
        DO I=1,N
                SUM1=SUM1+W(I)
        ENDDO 
        
        TB=SUM1/REAL(N)
        SUM1=0.

        DO I=1,N
                SUM1=SUM1+(W(I)-TB)**2
!               SUM1=SUM1+(W(I))**2
        ENDDO

        SX=SQRT((1/REAL(N))*SUM1)
        SX=SX/3.0
        Q33=(TB-SX)
        Q66=(TB+SX)

        RETURN
        END

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

      SUBROUTINE CHON_TINH

      INCLUDE 'PLOP.INC'
!       PARAMETER (MVAR1=8,NTOHOP=400, NNHOM=3, NMAX=200, MMAX=12)
!       PARAMETER (MVAR1=8,NTOHOP=400, NNHOM=3)
      character *70 fsl,fra,FTEMP
        INTEGER NGROUP, NVAR, IPT
        PARAMETER (NGROUP=NNHOM, NVAR=MVAR1) ! SO NHOM, SO BIEN, 
                                             ! SO MAU CUC DAI MOI NHOM
        REAL PTPL(NNHOM,MVAR1+1), DCX
      REAL PTR(NTOHOP,NNHOM,MVAR1+1), ACCURATE(NTOHOP) ! LUU CAC PHUONG
TRINH
                                                                                                  ! VA
                                                                                                  ! DO
                                                                                                  ! CHINH
                                                                                                  ! XAC

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

        REAL PP0(NTOHOP,NNHOM,NNHOM),P10(NTOHOP,NNHOM),P20(NTOHOP,NNHOM)
        COMMON /TABLE/ PP0,P10,P20
        COMMON /STORE_KQ/ PTR, ACCURATE, IPT
        COMMON /ACCU/ PTPL,DCX
        COMMON /FNAME/ fsl,fra,FTEMP
        INTEGER I1,I2,I3,I4,I5,I6,I7,I8,I9,I0
        COMMON /VARIDX/ I1,I2,I3,I4,I5,I6,I7,I8,I9,I10

      IPT=0   ! DEM SO PHUONG TRINH
        I1=0
        I2=0
        I3=0
        I4=0
        I5=0
        I6=0
        I7=0
        I8=0
        I9=0
        I10=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)

      INCLUDE 'PLOP.INC'
!       PARAMETER (MVAR1=8,NTOHOP=400, NNHOM=3, NMAX=200, MMAX=12)
!       PARAMETER (MVAR1=8,NTOHOP=400, NNHOM=3)
        REAL PTPL(NNHOM,MVAR1+1), DCX
      REAL PTR(NTOHOP,NNHOM,MVAR1+1), ACCURATE(NTOHOP) 
                            ! LUU CAC PHUONG TRINH
                                                ! VA DO CHINH XAC
        INTEGER I1,I2,I3,I4,I5,I6,I7,I8,I9,I10
        INTEGER IDXVAR(MVARMAX+1)

        REAL PP0(NTOHOP,NNHOM,NNHOM),P10(NTOHOP,NNHOM),P20(NTOHOP,NNHOM)
        REAL PP(NNHOM,NNHOM), P1(NNHOM),P2(NNHOM)   ! LUU KET QUA DANG
GIA
        COMMON /BANGDG/ PP,P1,P2

        COMMON /VARIDX/ I1,I2,I3,I4,I5,I6,I7,I8,I9,I10
        COMMON /STORE_KQ/ PTR, ACCURATE
        COMMON /TABLE/ PP0,P10,P20
        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
        IDXVAR(10)=I9 + 1
        IDXVAR(11)=I10 + 1
!       PRINT '(10I5)',(IDXVAR(I),I=1,9)
!       PAUSE
        PP0(IPT,:,:) = PP(:,:)
        P10(IPT,:  ) = P1(:  )
        P20(IPT,:  ) = P2(:  )

        DO I=1,NNHOM
           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
      INCLUDE 'PLOP.INC'
!       PARAMETER (MVAR1=8,NTOHOP=400, NNHOM=3, NMAX=200, MMAX=12)
!       PARAMETER (MVAR1=8,NTOHOP=400, NNHOM=3)
!       PARAMETER (NLOP=NNHOM)
      dimension N(NNHOM),CMEAN(MVAR1),XBAR(MVAR1*NNHOM),
        &          C((MVAR1+1)*NNHOM),D(MVAR1*MVAR1)
      DIMENSION P(NTOHOP*NNHOM),X(NTOHOP*NNHOM*MVAR1)
        integer   LG(NTOHOP*NNHOM)
      character *70 fsl,FTEMP
        DIMENSION ND(NNHOM),ND1(NNHOM),ND2(NNHOM),ND3(NNHOM)
        REAL PP(NNHOM,NNHOM), P1(NNHOM),P2(NNHOM)   ! LUU KET QUA DANG
GIA
        REAL PTPL(NNHOM,MVAR1+1), DCX

        COMMON /ACCU/ PTPL,DCX
      COMMON /FNAME/ fsl
        COMMON /DLMAU/ N
        COMMON /BANGDG/ PP,P1,P2
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 
!  190  write(3,14) L,P(J),LG(j)

   14           format(1x,i4,2x,F8.5,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
!######################################################