Real, Parameter :: Miss=-99.

      Real*8, Allocatable :: Xbar(:), Std(:) ! D(:), V(:)
      Real*8, Allocatable :: X(:,:), Y(:,:), Yest(:,:), YT(:,:), Xorg(:,:)
      Real, Allocatable :: Xnew(:), Fnew(:)
      ! Jacobi:
      Real, Allocatable :: RR(:,:)
      
      Real*8, Allocatable :: Fx(:,:), Ax(:,:)   ! X = Fx*A'x
      Real, Allocatable :: EOF(:,:), PCs(:,:)
      !For testing:
      Real, Allocatable :: FF(:,:), AA(:,:)
      Real, Allocatable :: X1(:,:), EVtor2(:,:)
      Real, Allocatable :: YdX(:,:)

      Real, Allocatable :: EVtor(:,:), Lamda(:,:)
      Real*8, Allocatable :: EValues(:)
      ! New:
      double precision, Allocatable :: Cor(:,:), Evect(:,:)
      double precision :: abserr

      Integer :: Nt, Nx
      Integer :: Mcount, Ncount, switch, Standard, Nmode, Nthre, PC_Stard
      Real :: Threshold
      Character (Len=250) :: Fname
      Real :: Factor

      Namelist /Model_Info/ Fname, switch, standard, PC_Stard, Nmode 

      Integer :: i,j,k

      Open (1,file="Namelist_PCA_New.txt",Status="old")
      Read (1, Model_Info)

      Nt=10; Nx=15
      Factor = Sqrt(Real(Nt))

      Allocate (X(Nt,Nx), Y(Nt,Nx), YT(Nx,Nt))   ! Origin data, Centered or standardized data and its transpose
      Allocate (Yest(Nt,Nx))
      Allocate (Xorg(Nt,Nx), Xnew(Nx),Fnew(Nmode))

      Open (1,file=fname)
      read (1,*)
      Do i=1, Nt
         read (1,*) (X(i,j),j=1,Nx)
      Enddo
      
      if (switch==1) then
         Mcount=Nt; Ncount=Nx
      else if (switch==0) then
         Mcount=Nx; Ncount=Nt
      endif
      !For testing:
      Allocate (FF(Mcount,Mcount), AA(Mcount,Mcount))
      Allocate (X1(Nt,Nx), EVtor2(Mcount,Mcount))
      Allocate (YdX(Nt,Nx))

      Allocate (Xbar(Mcount), Std(Mcount))                      ! TB, DLC
      Allocate (EValues(Mcount))                                ! Tri rieng va Vecto rieng (1D)
      Allocate (Fx(Ncount,Mcount), Ax(Mcount,Mcount))           ! EOF va PCA
      Allocate (EVtor(Mcount,Mcount))                ! Working arr va matran Vecto rieng
      Allocate (Lamda(Mcount,Mcount))
      ! Jacobi:
      Allocate (RR(Mcount,Mcount))
      ! Estimate of Fx and Ax:
      Allocate (EOF(Ncount,Mcount), PCs(Mcount,Mcount))           ! EOF va PCA

      Call PCA_2019 (X,Y,Nt,Nx,switch, Standard, PC_Stard, &
                     Xbar,Std,EValues,Fx,Ax,Nmode)

      K=Nmode

      Do i=1,Nt
         Do j=1,Nx
            Yest(i,j)=X(i,j)-Xbar(j)
            if (Standard==1) Yest(i,j)=Yest(i,j)/Std(j)
         Enddo
         print"(30F5.1)", (Y(i,j),Yest(i,j),j=1,Nx)
      Enddo
      pause
      Do j=1,Mcount
         lamda(j,j) = 1./EValues(j)            ! 1/sqrt(eig)
      Enddo
      EOF = Matmul(Matmul(Yest,Ax),lamda)/Factor
      !Do j=1,Mcount
      !   EOF(:,j)=EOF(:,j)/Factor
      !Enddo

!      Call Calc_EOF_PCA
      FF=MatMul(Transpose(Fx),Fx)
      AA=MatMul(Transpose(Ax),Ax)
      if (switch==0) then
         X1=MatMul(Fx,Transpose(Ax))*Factor
      else if (switch==1) then
         X1=MatMul(Fx,Transpose(Ax))*Factor
         !Do j=1,Nx
         !   Do i=1,Nt
         !      YdX(i,j)=Y(i,j)/X1(i,j)
         !   Enddo
         !Enddo
      endif
      ! Restore the origin values of X:
      
         Do j=1,Nx
            Do i=1,Nt
               if (standard==1) then
                  Xorg(i,j)=X1(i,j)*Std(j)+Xbar(j)
               else
                  Xorg(i,j)=X1(i,j)+Xbar(j)
               endif 
            Enddo
         Enddo
      
      !if (standard==1) then
      !   Do j=1, Mcount
         !   X1(:,j)=X1(:,j)*Factor
      !   Enddo
      !endif

      !EVtor2=MatMul(Transpose(EVtor),EVtor)

      print*,EValues
      print*,"Sum of Cova=",Sum(EValues)
      Call Write_test


      Contains
      !====================================================
      Include "/home/tanpv/Mylibs/Sub_PCA_2019.f90"

      !---------------------------------------------------
      Subroutine Write_test
      print*,"PCA"
      Do i=1,Nx
         print "(20F10.2)", (Ax(i,j),j=1,K) ! Mcount)
      Enddo
      pause
      print*,"EOF"
      Do i=1, Nt
         print "(20F10.2)", (Fx(i,j),j=1,K) ! Mcount)
      Enddo
      pause
      print*,"EOF Estimated:"
      Do i=1, Nt
         print "(20F10.2)", (EOF(i,j),j=1,K) ! Mcount)
      Enddo
      pause

      print*,"AA"
      Do i=1, K !Mcount
         print "(20F10.2)", (AA(i,j),j=1,K) ! Mcount)
      Enddo
      pause
      print*,"FF"
      Do i=1, K
         print "(20F10.2)", (FF(i,j),j=1,K)! Mcount)
      Enddo
      pause
      print*,"Y"
      Do i=1, Ncount
         print "(20F10.2)", (Y(i,j),j=1, Mcount)
      Enddo
      pause
      print*,"X1"
      Do i=1, Ncount
         print "(20F10.2)", (X1(i,j),j=1, Mcount)
      Enddo
      pause
      print*,"Xorg"
      Do i=1, Ncount
         print "(20F10.2)", (X(i,j),j=1, Mcount)
      Enddo
      pause
      print*,"X-estimated"
      Do i=1, Ncount
         print "(20F10.2)", (Xorg(i,j),j=1, Mcount)
      Enddo
      pause
      !print*,"EVtor2"
      !Do i=1, Mcount
      !   print "(20F10.2)", (EVtor2(i,j),j=1, Mcount)
      !Enddo
      End Subroutine Write_test
      !---------------------------------------------------

      END