$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
!######################################################