! Chuong trinh tinh PCs, EOFs theo cach moi (2019-10-10)
      ! Input:
      !   X(Nt,Nx): Mang so lieu kich thuoc Nt (Time), Nx (Space)
      !             X(Nt,Nx) = EOFs(Nt,K).PCs'(K,Nx)
      !   switch: Tham so de quyet dinh su dung cach tinh ma tran tuong quan:
      !   if (switch==1) then
      !       Mcount=Nt; Ncount=Nx
      !       Khi do: X(Nx,Nt) = PCs(Nx,K).EOFs'(K,Nt)
      !   else if (switch==0) then
      !       Mcount=Nx; Ncount=Nt
      !       Khi do: X(Nt,Nx) = EOFs(Nt,K).PCs'(K,Nx)
      !   endif
      !   Standard: =1: su dung ma tran Correlation, =0: ma tran Covariance
      !   PC_Stard: =1: chuan hoa cac PCs (Cac PCs truc giao chuan hoa)
      !   Nmode: So mode duoc giu lai (Nmode=K)
      ! Output:
      !   Xbar(Nx), Std(Nx): Trung binh va do lech chuan cua X tinh voi Nt la dung luong mau
      !   Y(Nt,Nx)=X(Nt,Nx)-Xbar(Nx) neu PC_Stard=1
      !   Y(Nt,Nx)=(X(Nt,Nx)-Xbar(Nx))/Std(Nx) neu PC_Stard=0
      !   EValues(Nmode): Cac gia tri rieng
      !   EOF(Nt,Nmode), PCs(Nx,Nmode): Cac ham thoi gian va khong gian cua X(Nt,Nx)
      !
      Subroutine PCA_2019 (X,Y,Nt,Nx,switch, Standard, PC_Stard, &
                           XBAR,STD,EValues,EOF,PCs,Nmode)
        Integer :: Nt,Nx,Ncount,Mcount,switch, Standard, PC_Stard,Nmode
        Integer :: i,j,K
        Real :: SQ, Factor
        Real*8 :: X(Nt,Nx)
        Real*8, Allocatable :: Y(:,:), Xbar(:), Std(:)
        Real*8, Allocatable :: EOF(:,:), PCs(:,:), Fx(:,:), Ax(:,:)
        Real*8, Allocatable :: RR(:,:), EVtor(:,:)
        Real*8, Allocatable :: EValues(:)

        if (switch==1) then
           Mcount=Nt; Ncount=Nx
        else if (switch==0) then
           Mcount=Nx; Ncount=Nt
        endif
        Factor = Sqrt(Real(Nt))
        if (.not.Allocated (Y))    Allocate (Y(Nt,Nx))
        if (.not.Allocated (Xbar)) Allocate (Xbar(Nx))
        if (.not.Allocated (Std))  Allocate (Std(Nx))
        if (.not.Allocated (RR))   Allocate (RR(Mcount,Mcount))   ! TB, DLC, ma tran Covar hoac Corre
        if (.not.Allocated (EValues)) Allocate (EValues(Mcount))
        if (.not.Allocated (EVtor))   Allocate (EVtor(Mcount,Mcount))          ! Tri rieng va Vecto rieng (1D)
        if (.not.Allocated (Fx))      Allocate (Fx(Ncount,Mcount)) 
        if (.not.Allocated (Ax))      Allocate (Ax(Mcount,Mcount))           ! EOF va PCA
        if (.not.Allocated (EOF))     Allocate (EOF(Nt,Nmode)) 
        if (.not.Allocated (PCs))     Allocate (PCs(Nx,Nmode))
        !---------------------------------------------------------------------------------

        Call Covar (X,Y,Nt,Nx,Mcount,RR,XBAR,STD,switch,standard)
        Call Eigen (RR,Mcount,EValues,EVtor)
        Call eigsrt (EValues,EVtor,Mcount)
         K = Nmode
         DO j=1,K
            SQ=Sqrt(EValues(j))
            DO i=1,Mcount
               Ax(I,J)=EVtor(i,j)*SQ
            ENDDO
         ENDDO
         if (PC_Stard==1) then
            Call Calc_EOF_PCA_Standard (Y,Nt,Nx,Ax,EVtor,EValues,Ncount,Mcount,EOF,PCs,Nmode,switch)
         else if (PC_Stard==0) then
            Call Calc_EOF_PCA (Y,Nt,Nx,Ax,EValues,Ncount,Mcount,EOF,PCs,Nmode,switch)
         endif
      End Subroutine PCA_2019
      !----------------------------------------------------------------------------------!
      Subroutine Calc_EOF_PCA_Standard (Y,Nt,Nx,Ax,EVtor,EValues,Ncount,Mcount,EOF,PCs,Nmode,switch)
      Integer :: i,j,K,Nt,Nx,Ncount,Mcount,switch,Nmode
      Real*8 :: Modul_Fx, Factor
      Real*8, Allocatable :: Lamda(:,:), YT(:,:), Fx(:,:)
      Real*8 :: Y(Nt,Nx) !, Ax(Mcount,Mcount), EVtor(Mcount,Mcount)
      Real*8 :: Ax(Mcount,Nmode), EVtor(Mcount,Nmode), EValues(Nmode)
      Real*8 :: EOF(Nt,Nmode), PCs(Nx,Nmode)
      !Real*8 :: EValues(Mcount)

        !Allocate (Lamda(Mcount,Mcount),YT(Nx,Nt),Fx(Ncount,Mcount))
        K = Nmode
        Allocate (Lamda(K,K),YT(Nx,Nt),Fx(Ncount,K))
        Factor = Sqrt(Real(Nt))
        if (switch==0) then
           Ax(1:Mcount,1:K)=EVtor(1:Mcount,1:K)
           Fx(1:Ncount,1:K) = MatMul(Y(1:Ncount,1:Mcount),Ax(1:Mcount,1:K))
           Do j=1,K
              Fx(:,j)=Fx(:,j)/Factor               ! Voi PCs chuan hoa (Ax=EVtor)
           Enddo
           EOF(1:Ncount,1:K) = Fx(1:Ncount,1:K)
           PCs(1:Mcount,1:K) = Ax(1:Mcount,1:K)
        else if (switch==1) then
           Do j=1, K
              Lamda(j,j)=Sqrt(EValues(j))
           Enddo
           YT = Transpose(Y)
           Fx(1:Ncount,1:K) = MatMul(YT(1:Ncount,1:Mcount),EVtor(1:Mcount,1:K))
           Do j=1,K
              Modul_Fx=0
              Do i=1, Ncount
                 Modul_Fx = Modul_Fx + Fx(i,j)*Fx(i,j)
              Enddo
              Modul_Fx = Sqrt(Modul_Fx)
              Fx(:,j)=Fx(:,j)/Modul_Fx
           Enddo
           EOF(1:Mcount,1:K) = Ax(1:Mcount,1:K)
           PCs(1:Ncount,1:K) = Fx(1:Ncount,1:K)
        endif
        Deallocate (Lamda)
      End Subroutine Calc_EOF_PCA_Standard
      !---------------------------------------------------------------------------------
      Subroutine Calc_EOF_PCA (Y,Nt,Nx,Ax,EValues,Ncount,Mcount,EOF,PCs,Nmode,switch)
      Integer :: i,j,K,Nt,Nx,Ncount,Mcount,switch,Nmode
      Real*8 :: Modul_Fx, Factor, SQE
      Real*8, Allocatable :: Lamda(:,:), YT(:,:), Fx(:,:)
      Real*8 :: Y(Nt,Nx) !, Ax(Mcount,Mcount)
      Real*8 :: Ax(Mcount,Nmode), EValues(Nmode)
      Real*8 :: EOF(Nt,Nmode), PCs(Nx,Nmode)
      !Real*8 :: EValues(Mcount)
      
        !Allocate (Lamda(Mcount,Mcount),YT(Nx,Nt),Fx(Ncount,Mcount))
        K = Nmode
        Allocate (Lamda(K,K),YT(Nx,Nt),Fx(Ncount,K))
        Factor = Sqrt(Real(Nt))
        if (switch==0) then
           Fx(1:Ncount,1:K) = MatMul(Y(1:Ncount,1:Mcount),Ax(1:Mcount,1:K))
           Do j=1,K
              Fx(:,j)=Fx(:,j)/(EValues(j))/Factor
           Enddo
           EOF(1:Ncount,1:K) = Fx(1:Ncount,1:K)
           PCs(1:Mcount,1:K) = Ax(1:Mcount,1:K)
        else if (switch==1) then
           Do j=1, K
              SQE = sqrt(EValues(j))
              Ax(:,j)=Ax(:,j)/SQE  ! sqrt(EValues(j))
              Lamda(j,j)=SQE       ! Sqrt(EValues(j))
           Enddo
           YT = Transpose(Y)
           Fx(1:Ncount,1:K) = MatMul(YT(1:Ncount,1:Mcount),Ax(1:Mcount,1:K))
           Do j=1,K
              Modul_Fx=0
              Do i=1, Ncount
                 Modul_Fx = Modul_Fx + Fx(i,j)*Fx(i,j)
              Enddo
              Modul_Fx = Sqrt(Modul_Fx)
              Fx(:,j)=Fx(:,j)/Modul_Fx
           Enddo
           Fx=MatMul(Fx(1:Ncount,1:K),Lamda(1:K,1:K))
           EOF(1:Mcount,1:K) = Ax(1:Mcount,1:K)
           PCs(1:Ncount,1:K) = Fx(1:Ncount,1:K)
        endif
        Deallocate (Lamda)
        print*,"Calculating EOF & PCs ... Done !"
      End Subroutine Calc_EOF_PCA
      !---------------------------------------------------------------------------------
      SUBROUTINE Covar(X,Y,Nt,Nx,Mcount,RR,XBAR,STD,switch,standard)
      Integer :: Nt, Nx, Mcount, N1, L, J, K, IO, I,switch,standard
      Real*8 :: Sx, Xtb
      Real*8 :: X(Nt,Nx), Y(Nt,Nx), Xbar(Nx), Std(Nx), RR(Mcount,Mcount)
      !Real*8 :: YT(Nx,Nt)
      !---------------------------------------------------
      Do j=1, Nx
         Sx=0; Xtb=0; N1=0
         Do i=1, Nt
            if (X(i,j)/=Miss) then
               n1=n1+1; Xtb=Xtb+X(i,j); Sx=Sx+X(i,j)*X(i,j)
            Endif
         Enddo
         if (n1>nt/2) then
            XBAR(j)=Xtb/real(n1); STD(j)=Sqrt(Sx/real(n1)-XBAR(j)*XBAR(j))
         endif
         Do i=1, Nt
            Y(i,j) = (X(i,j)-XBAR(j))
            if (standard==1) Y(i,j) = Y(i,j)/STD(j)
         Enddo
      Enddo
      !YT=Transpose(Y)
      if (switch==1) then
         print"(30F6.2)",Y(Nt-30:Nt,Nx-30:Nx)
         print*,"Dim(RR)=",Nt,Mcount; pause
         !RR = MatMul(Y,YT)
         RR = MatMul(Y,Transpose(Y))
         !Call Y_YT(Y,YT,RR,Nt,Nx)
         print"(30F6.2)",RR(1:30,1:30); pause
      else if (switch==0) then
         !RR = MatMul(YT,Y) 
         RR = MatMul(Transpose(Y),Y)
         !Call YT_Y(YT,Y,RR,Nx,Nt)
      endif
      RR(:,:)=RR(:,:)/real(Nt)
      print*,"Calculating Covariance ... Done !"
      print"(30F6.2)",RR(1:30,1:30); pause
      End SUBROUTINE Covar
      !---------------------------------------------------------------------------------
      Subroutine Y_YT(Y,YT,RR,Nt,Nx)
         Integer :: i,j,k,Nt,Nx
         Real*8 :: Y(Nt,Nx),YT(Nx,Nt),RR(Nt,Nt) 

         Do j=1,Nt
            Do k=1,Nt
               RR(j,k)=0.0
               Do i=1,Nx
                  RR(j,k)=RR(j,k)+Y(j,i)*YT(i,k)
!print*,RR(j,k)
               Enddo
            Enddo
         Enddo
      End Subroutine Y_YT
      !---------------------------------------------------------------------------------
      Subroutine YT_Y(YT,Y,RR,Nx,Nt)
         Integer :: i,j,k,Nt,Nx
         Real*8 :: YT(Nx,Nt),Y(Nt,Nx),RR(Nx,Nx)       

         Do j=1,Nx
            Do k=1,Nx
               RR(j,k)=0.0
               Do i=1,Nt
                  RR(j,k)=RR(j,k)+YT(j,i)*Y(i,k)
               Enddo
            Enddo
         Enddo
      End Subroutine YT_Y
      !---------------------------------------------------------------------------------
      SUBROUTINE eigsrt (EValues,EVtor,Mcount)
      INTEGER :: n,np
      Real*8, Allocatable :: D(:), V(:,:)
      INTEGER :: i,j,k,Mcount
      REAL*8 :: p
      Real*8 :: EVtor(Mcount,Mcount)
      Real*8 :: EValues(Mcount)

      Allocate (D(Mcount),V(Mcount,Mcount))
      N=Mcount
      D(:)=EValues(:); V(:,:)=EVtor(:,:)

      do  i=1,n-1
        k=i
        p=d(i)
        do  j=i+1,n
          if(d(j).ge.p)then
            k=j
            p=d(j)
          endif
        enddo
        if(k.ne.i)then
          d(k)=d(i)
          d(i)=p
          do  j=1,n
            p=v(j,i)
            v(j,i)=v(j,k)
            v(j,k)=p
          enddo
        endif
      enddo
      EValues(:)=D(:)
      EVtor(:,:)=V(:,:)
      Deallocate (D, V)
      END SUBROUTINE eigsrt
      !---------------------------------------------------------------------------------
      SUBROUTINE Eigen (RR,Mcount,EValues,EVtor) 
      INTEGER :: n,np,nrot 
      INTEGER :: i,ip,iq,j
      REAL*8 :: c,g,h,s,sm,t,tau,theta,tresh 
      Integer :: Mcount
      Real*8 :: RR(Mcount,Mcount), EVtor(Mcount,Mcount)
      Real*8 :: EValues(Mcount)

      Real*8, Allocatable :: A(:,:), D(:), V(:,:), B(:), Z(:)

      Allocate (A(Mcount,Mcount),D(Mcount),V(Mcount,Mcount),B(Mcount),Z(Mcount))
      A(:,:) = RR(:,:); N=Mcount
      do  ip=1,n
        do  iq=1,n
          v(ip,iq)=0.
        enddo
        v(ip,ip)=1.
      enddo
      do  ip=1,n
        b(ip)=a(ip,ip)
        d(ip)=b(ip)
        z(ip)=0.
      enddo
      nrot=0
      do  i=1,50
        sm=0.
        do  ip=1,n-1
          do iq=ip+1,n
            sm=sm+abs(a(ip,iq))
          enddo
        enddo
        if(sm.eq.0.) then
          EValues(:)=D(:)
          EVtor(:,:)=V(:,:)
          Deallocate (A, D, V, B, Z)
          return
        endif
        if(i.lt.4)then
          tresh=0.2*sm/n**2
        else
          tresh=0.
        endif
        do  ip=1,n-1
          do  iq=ip+1,n
            g=100.*abs(a(ip,iq))
            if((i.gt.4).and.(abs(d(ip))+g.eq.abs(d(ip))).and. &
               (abs(d(iq))+g.eq.abs(d(iq))))then
              a(ip,iq)=0.
            else if(abs(a(ip,iq)).gt.tresh)then
              h=d(iq)-d(ip)
              if(abs(h)+g.eq.abs(h))then
                t=a(ip,iq)/h
              else
                theta=0.5*h/a(ip,iq)
                t=1./(abs(theta)+sqrt(1.+theta**2))
                if(theta.lt.0.)t=-t
              endif
              c=1./sqrt(1+t**2)
              s=t*c
              tau=s/(1.+c)
              h=t*a(ip,iq)
              z(ip)=z(ip)-h
              z(iq)=z(iq)+h
              d(ip)=d(ip)-h
              d(iq)=d(iq)+h
              a(ip,iq)=0.
              do  j=1,ip-1
                g=a(j,ip)
                h=a(j,iq)
                a(j,ip)=g-s*(h+g*tau)
                a(j,iq)=h+s*(g-h*tau)
              enddo
              do  j=ip+1,iq-1
                g=a(ip,j)
                h=a(j,iq)
                a(ip,j)=g-s*(h+g*tau)
                a(j,iq)=h+s*(g-h*tau)
              enddo
              do  j=iq+1,n
                g=a(ip,j)
                h=a(iq,j)
                a(ip,j)=g-s*(h+g*tau)
                a(iq,j)=h+s*(g-h*tau)
              enddo
              do  j=1,n
                g=v(j,ip)
                h=v(j,iq)
                v(j,ip)=g-s*(h+g*tau)
                v(j,iq)=h+s*(g-h*tau)
              enddo
              nrot=nrot+1
            endif
          enddo
        enddo
        do  ip=1,n
          b(ip)=b(ip)+z(ip)
          d(ip)=b(ip)
          z(ip)=0.
        enddo
      enddo
      Deallocate (A, D, V, B, Z)
      !pause 'too many iterations in jacobi'
      END SUBROUTINE Eigen 
      !----------------------------------------------------------------------------------!