! Chuong trinh tinh PCs, EOFs theo cach moi (2019-10-10) Use NetCDF Real, Parameter :: Miss=-99. Real, Allocatable :: Xbar(:), Std(:), Rx(:), D(:), V(:) Real, Allocatable :: X(:,:), Y(:,:), YT(:,:) ! X(nt,nx): SL ban dau, Y(nt,nx): SL quy tam hoac chuan hoa, YT-Chuyen vi cua Y Real, Allocatable :: Fx(:,:), Ax(:,:) ! X = Fx*A'x ! Fx, Ax la cac mang tinh trung gian Real, Allocatable :: EOF(:,:), PCs(:,:) ! EOF(nt,K), PCs(nx,K): Cac ham thoi gian va ham khong gian Real, Allocatable :: EVtor(:,:), Lamda(:,:) Real, Allocatable :: EValues(:), ECumulate(:) !--------------------------------------------------------------------------------- !For PCA's calculations: Real, Allocatable :: RR(:,:), A(:), R(:) Integer :: Nt, Nx ! Kich thuoc thuc cua mang X Integer :: Mcount, Ncount, switch, Standard, Nth0, Nthre, K Real :: Threshold Character (Len=250) :: Fname, Fname1 Real :: Factor Namelist /Model_Info/ Fname, switch, standard, Nmode Integer :: i,j,L ! Mo ta thong tin SL dau vao: !--------------------------------------------------------------------------------- Real, allocatable :: xlon(:), xlat(:), xtime(:) Integer :: Nlon, Nlat, Ntime Logical :: Status Real, allocatable :: PDSI(:,:,:), SPI12(:,:,:), SPEI12(:,:,:), Mask(:,:) Real, allocatable :: Dr_IDX(:,:,:) Integer :: Nst, iidx Integer, Parameter :: NDr_idx=3 ! PDSI, SPEI12, SPI12 Character (Len=5), Parameter :: IDX_nm(NDr_idx)= & (/"pdsi","spi12","spe12"/) !(/"PDSI","SPI12","SPE12"/) !--------------------------------------------------------------------------------- !For Writing EOFs, PCs: Real, allocatable :: PCA(:,:,:), xmode(:), Xtmp(:), EOF1(:,:) !--------------------------------------------------------------------------------- Open (1,file="Namelist_PCA_New.txt",Status="old") Read (1, Model_Info) Close (1) !--------------------------------------------------------------------------------- ! Reading data: Call Read_Drough_IDX (Fname) Nt = Ntime; Nx = Nst ! Danh cho qua trinh tinh EOFs, PCs: !--------------------------------------------------------------------------------- if (switch==1) then Mcount=Nt; Ncount=Nx else if (switch==0) then Mcount=Nx; Ncount=Nt endif Factor = Sqrt(Real(Nt)) Allocate (X(Nt,Nx), Y(Nt,Nx), YT(Nx,Nt)) ! Loop over all drought indices: !Fname="pca_drou_idx.nc"; Fname1="eof_timeseries_drou_idx.nc" !Do iidx=1, NDr_idx-1 Allocate (Xbar(Mcount), Std(Mcount), Rx(Mcount*Mcount)) ! TB, DLC, ma tran Covar hoac Corre Allocate (EValues(Mcount), V(Mcount*Mcount)) ! Tri rieng va Vecto rieng (1D) Allocate (Fx(Ncount,Mcount), Ax(Mcount,Mcount)) ! EOF va PCA Allocate (D(Mcount), EVtor(Mcount,Mcount)) ! Working arr va matran Vecto rieng Allocate (Lamda(Mcount,Mcount)) Allocate (EOF(Nt,Nmode), PCs(Nx,Nmode), ECumulate(Mcount)) !--------------------------------------------------------------------------------- !For PCA's calculations: Allocate (RR(Mcount,Mcount), A(Mcount*Mcount), R(Mcount*Mcount)) !--------------------------------------------------------------------------------- !For Writing EOFs, PCs: Allocate (PCA(Nlon,Nlat,Nmode), xmode(Nmode), Xtmp(Nmode), EOF1(Nmode,Nt)) !--------------------------------------------------------------------------------- ! Loop over all drought indices: Fname="pca_drou_idx.nc"; Fname1="eof_timeseries_drou_idx.nc" Do iidx=1,NDr_idx+1 if (iidx>NDr_idx) exit X(1:Nt,1:Nx) = Dr_IDX(1:Nt,1:Nx,iidx) ! Tinh cac EOFs, PCs Call Covar(Nt,Nx, Miss) Call EIGEN(Mcount,MV) Call Calc_Thresholds (Nmode, CON) Call TRACE(Mcount,Nmode,K) Call EigenValues_and_EigenVectors CALL LOAD(Mcount,K) Call Calc_AX Call Calc_EOF_PCA !------------------------------------------------------------------------------ ! Write results to files: ! Fname="pca_drou_idx.nc"; Fname1="eof_timeseries_drou_idx.nc" Call Write_EOF_PCA (Fname,Fname1,iidx) print*, "===========> Complete for ",IDX_nm(iidx) Enddo !Deallocate(Xbar,Std,Rx,EValues,V,Fx,Ax,D,EVtor,Lamda,EOF,PCs,ECumulate) !Enddo print*,"Program is normal terminating ..." CONTAINS !--------------------------------------------------------------------------------- Subroutine Read_Drough_IDX (Fname) Integer :: ist, ilon, ilat, it, Nt, imon, imon1, iyr Integer :: Nyear, Nmon, nlat1, nlat0 Character (Len=*) :: Fname Real :: tmp, xmean(12), ic(12) Real, Allocatable :: X0(:), X1(:,:,:), X2(:,:,:), X3(:,:,:) Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, PDSI, "ori_pdsi") Deallocate (xlon, xlat, xtime) Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, SPEI12, "SPEI_12") Deallocate (xlon, xlat, xtime) Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, SPI12, "SPI_12") Deallocate (xlon, xlat) Call rd_ncdf_2d (trim(Fname), xlon,xlat, nlon,nlat,Mask,"land") Nyear=30 ! 1976-2005 Nmon =12 Ntime = Nt - Nmon Nlat1 = nlat Nlat0 = 0 !15 Nlat = nlat-Nlat0 ! Tru vung Australia (1-15 ~ 15S-11S) Allocate (X0(nlat)); X0(1:nlat)=xlat(Nlat0+1:Nlat1) Deallocate (xlat); Allocate (xlat(nlat)) xlat(:) = X0(:) Deallocate (X0) Allocate (X1(nlon, nlat, Ntime),X2(nlon, nlat, Ntime),X3(nlon, nlat, Ntime)) X1(:,1:nlat,1:Ntime)=PDSI(:,Nlat0+1:Nlat1,Nmon+1:Nt) X2(:,1:nlat,1:Ntime)=SPEI12(:,Nlat0+1:Nlat1,Nmon+1:Nt) X3(:,1:nlat,1:Ntime)=SPI12(:,Nlat0+1:Nlat1,Nmon+1:Nt) Deallocate (PDSI,SPEI12,SPI12) Allocate (PDSI(nlon,nlat,Ntime),SPEI12(nlon,nlat,Ntime),SPI12(nlon,nlat,Ntime)) PDSI=X1; SPEI12=X2; SPI12=X3 Deallocate (X1, X2, X3) ist=0 Do ilat=1, Nlat Do ilon=1, Nlon if (Mask(ilon,ilat)>0) ist=ist+1 Enddo Enddo Nst=ist Allocate (Dr_IDX(Ntime,Nst,NDr_idx)) ist=0 Do ilat=1, Nlat Do ilon=1, Nlon if (Mask(ilon,ilat)>0) then ist=ist+1 Dr_IDX(1:Ntime,ist,1)=PDSI (ilon,ilat,1:Ntime) Dr_IDX(1:Ntime,ist,2)=SPI12 (ilon,ilat,1:Ntime) Dr_IDX(1:Ntime,ist,3)=SPEI12(ilon,ilat,1:Ntime) endif Enddo Enddo Deallocate (PDSI,SPEI12,SPI12) print*,"Reading data ... Done !" End Subroutine Read_Drough_IDX !--------------------------------------------------------------------------------- Subroutine Write_EOF_PCA (F_pca, F_eof, iidx) Character (Len=*) :: F_eof, F_pca Character (Len=255) :: var_ln, t_unit Integer :: i,j,k,iidx !Real, Allocatable :: PCA(:,:,:), xmode(:), EOF1(:,:), X1(:) !------------------------------------------------------------------------- PCA(:,:,:)=Miss; k=0 Do j=1, Nlat Do i=1, Nlon if (Mask(i,j)>0) then k=k+1; PCA(i,j,1:Nmode) = PCs(k,1:Nmode) endif Enddo Enddo Do i=1, Nmode; xmode(i)=i; Enddo ! PCA's modes var_ln="Spatial distribution of PCAs of "//trim(IDX_nm(iidx)) t_unit="months since 1976-01-01 " if (iidx==1) then Call wr_ncdf_3d (F_pca,xlon,xlat,xmode,nlon,nlat,Nmode, PCA, & trim(IDX_nm(iidx)), trim(t_unit),trim(var_ln), & "-", "lon", "lat","mode", Miss) else Call wr_ncdf_3d_apd (F_pca,PCA,trim(IDX_nm(iidx)),nlon,nlat, & Nmode, trim(var_ln), "-", "lon", "lat", "mode", Miss) endif !------------------------------------------------------------------------- Xtmp(1:Nmode) = EValues(1:Nmode) print*,"Eigenvalues:",Xtmp Call wr_ncdf_1d_apd (F_pca, Xtmp, trim(IDX_nm(iidx))//"_eig", Nmode, & "Eigen values of "//trim(IDX_nm(iidx)), "-", "mode", Miss) !X(1:nmodes) = variance(1:nmodes) !Call wr_ncdf_1d_apd (F_pca, X, trim(IDX_nm(iidx))//"_var",nmodes, & ! "Variances of "//trim(IDX_nm(iidx)), "-", "mode", Miss ) Xtmp(1:Nmode) = ECumulate(1:Nmode) print*,ECumulate(1:Nmode) print*,Xtmp;pause Call wr_ncdf_1d_apd (F_pca, Xtmp, trim(IDX_nm(iidx))//"_cum",Nmode, & "Cumuluted Eigenvalues of "//trim(IDX_nm(iidx)), & "-", "mode", Miss ) !------------------------------------------------------------------------- Do i=1, Nmode; EOF1(i,1:Nt)=EOF(1:Nt,i); Enddo var_ln="Timeseries of eofs of "//trim(IDX_nm(iidx)) if (iidx==1) then Call wr_ncdf_2d (F_eof,xmode,xtime, Nmode,Nt, EOF1, "eof_"//trim(IDX_nm(iidx)), & trim(var_ln), "-", "mode", "mon", Miss) else Call wr_ncdf_2d_apd (F_eof,EOF1,"eof_"//trim(IDX_nm(iidx)), Nmode,Nt, & trim(var_ln), "-", "mode", "mon", Miss) endif print*,"Write_EOF_PCA: Done !!!" End Subroutine Write_EOF_PCA !--------------------------------------------------------------------------------- !--------------------------------------------------------------------------------- Subroutine Calc_AX Integer :: i,j DO I=1,Mcount DO J=1,K L=Mcount*(J-1)+I D(J)=V(L) Ax(I,J)=V(L) ! Ma tran PCs, moi cot 1 PC ENDDO ENDDO print*,"Calculating Ax ... Done !" End Subroutine Calc_AX !--------------------------------------------------------------------------------- Subroutine EigenValues_and_EigenVectors Integer :: i, L DO I=1,K L=I+(I*I-I)/2 EValues(I)=Rx(L) ENDDO L=0 DO J=1,K DO I=1,Mcount L=L+1 D(I)=V(L) EVtor(i,j)=V(L) ! Ma tran Vector rieng, moi cot 1 vector ENDDO ENDDO print*,"Calculating EigenValues_and_EigenVectors ... Done !" End Subroutine EigenValues_and_EigenVectors !--------------------------------------------------------------------------------- Subroutine Calc_Thresholds (Nmode, CON) ! Input: Number of Eigenvalues (Nth0) and Threshold value (Threshold) of Eigenvalues ! Output: CON satisfying N-eigenvalue>=Nth0 OR accumulated Eigenvalues>Threshold Real :: tmp, tmp1, CON Integer :: i, j, L, Nthre, Nmode tmp = 0. do I=1,Mcount L=I+(I*I-I)/2 EValues(I) = Rx(L) tmp = tmp+EValues(I) ! Tong cac tri rieng enddo Do i=1,Mcount ECumulate(i)=0 Do j=1,I ECumulate(i)=ECumulate(i) + EValues(j)/tmp Enddo !Enddo !Do i=1,Mcount ! ECumulate(i)=ECumulate(i)/tmp Enddo print*,"Tong tri rieng=",tmp print*,"====>>",ECumulate(1:Nmode) tmp1 = 0; Nthre = 0 ! De so voi Tong cac tri rieng Do i=1,Mcount tmp1 = tmp1+EValues(I); Nthre=Nthre+1 if (Nthre >= Nmode) exit !if (tmp1/tmp > Threshold .OR. Nthre >= Nth0) exit ! Vuot nguong nhung khong qua Nth0 nhan to Enddo CON = EValues(Nthre) ! Xac dinh lai nguong gia tri cua Tri rieng duoc giu lai print*,"Evalues=",EValues(1:Nmode) print*,"CON=",CON WRITE(*,*)' TRI RIENG 1,2,3,4,5=',Rx(1),Rx(3),Rx(6),Rx(10),Rx(15) write(*,'(A,F10.3)') "Gia tri nguong CON: ", CON write(*,'(A,F10.3)') "So phan tram luong thong tin duoc giu lai: ",tmp1/tmp*100. print*,"Calculating Threshold ... Done !" End Subroutine Calc_Thresholds !--------------------------------------------------------------------------------- Subroutine Calc_EOF_PCA Integer :: i,j Real :: Modul_Fx if (switch==0) then Fx(1:Nt,1:K) = MatMul(Y(1:Nt,1:Nx),Ax(1:Nx,1:K)) Do j=1,K Fx(:,j)=Fx(:,j)/(EValues(j))/Factor Enddo EOF(1:Nt,1:K) = Fx(1:Nt,1:K) PCs(1:Nx,1:K) = Ax(1:Nx,1:K) else if (switch==1) then Do j=1, K Ax(:,j)=Ax(:,j)/sqrt(EValues(j)) Lamda(j,j)=Sqrt(EValues(j)) Enddo YT = Transpose(Y) Fx(1:Nx,1:K) = MatMul(YT(1:Nx,1:Nt),Ax(1:Nt,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:Nt,1:K) = Ax(1:Nt,1:K) PCs(1:Nx,1:K) = Fx(1:Ncount,1:K) endif print*,"Calculating EOF & PCs ... Done !" End Subroutine Calc_EOF_PCA !--------------------------------------------------------------------------------- SUBROUTINE Covar(Nt,Nx, Miss) Integer :: Nt, Nx, N1, L, J, K, IO, I Real :: Miss, Sx, Xtb !--------------------------------------------------- Where (IsNaN(X)); X=Miss; Endwhere 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)) if (IsNaN(Xbar(j))) print*,"XBAR NaN",j endif Do i=1, Nt if (X(i,j)==Miss) then X(i,j)=XBAR(j) ! Fill-up Miss by mean values endif Y(i,j) = (X(i,j)-XBAR(j)) if (standard==1) Y(i,j) = Y(i,j)/STD(j) Enddo Enddo if (switch==1) then RR = MatMul(Y,Transpose(Y)) else if (switch==0) then RR = MatMul(Transpose(Y),Y) endif RR(:,:)=RR(:,:)/real(Nt) ! Store R in 1D-Array of Rx to be returned Do I=1,Mcount DO J=1,Mcount if ((I-J) < 0) then L=I+(J*J-J)/2 else L=J+(I*I-I)/2 endif RX(L)=RR(i,j) ENDDO ENDDO print*,"Calculating Covariance ... Done !" End SUBROUTINE Covar !--------------------------------------------------------------------------------- SUBROUTINE EIGEN(N,MV) ! Input: Covariance matrix of Rx(:) with real size is M*M elements ! switch: =1 neu tinh tuong quan theo cac lat cat thoi gian (dao chieu) ! Output: Rx with Eigenvalues are stored on diagonal ! V(:) matrix with real size is M*M elements stored Eigenvectors in column by column REAL ANORM,ANRMX,THR,X,Y,SINX,SINX2,COSX,COSX2,SINCS,XRANGE Integer :: N, I, J, M, icount !WRITE(*,*)' BAT DAU EIGEN, N=', N icount=0 A(:) = Rx(:) 5 XRANGE=1.E-6 IF(MV-1) 10,25,10 10 IQ=-N DO 20 J=1,N IQ=IQ+N DO 20 I=1,N IJ=IQ+I R(IJ)=0. IF(I-J) 20,15,20 15 R(IJ)=1. 20 CONTINUE 25 ANORM=0. DO 35 I=1,N DO 35 J=1,N IF(I-J) 30,35,30 30 IA=I+(J*J-J)/2 ANORM=ANORM+A(IA)*A(IA) 35 CONTINUE IF(ANORM)165,165,40 40 ANORM=1.414*SQRT(ANORM) ANRMX=ANORM*XRANGE/FLOAT(N) IND=0 THR=ANORM 45 THR=THR/FLOAT(N) !icount=icount+1; print*,"Lan lap: ",icount 50 L=1 55 M=L+1 60 MQ=(M*M-M)/2 LQ=(L*L-L)/2 LM=L+MQ 62 IF(ABS(A(LM))-THR)130,65,65 65 IND=1 LL=L+LQ MM=M+MQ X=.5*(A(LL)-A(MM)) 68 Y=-A(LM)/SQRT(A(LM)*A(LM)+X*X) IF(X) 70,75,75 70 Y=-Y 75 SINX=Y/SQRT(2.*(1.+(SQRT(1.-Y*Y)))) SINX2=SINX*SINX 78 COSX=SQRT(1.-SINX2) COSX2=COSX*COSX SINCS=SINX*COSX ILQ=N*(L-1) IMQ=N*(M-1) DO 125 I=1,N IQ=(I*I-I)/2 IF(I-L) 80,115,80 80 IF(I-M) 85,115,90 85 IM=I+MQ GOTO 95 90 IM=M+IQ 95 IF(I-L) 100,105,105 100 IL=I+LQ GOTO 110 105 IL=L+IQ 110 X=A(IL)*COSX-A(IM)*SINX A(IM)=A(IL)*SINX+A(IM)*COSX A(IL)=X 115 IF(MV-1) 120,125,120 120 ILR=ILQ+I IMR=IMQ+I X=R(ILR)*COSX-R(IMR)*SINX R(IMR)=R(ILR)*SINX+R(IMR)*COSX R(ILR)=X 125 CONTINUE X=2.*A(LM)*SINCS Y=A(LL)*COSX2+A(MM)*SINX2-X X=A(LL)*SINX2+A(MM)*COSX2+X A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) A(LL)=Y A(MM)=X 130 IF(M-N) 135,140,135 135 M=M+1 GOTO 60 140 IF(L-(N-1)) 145,150,145 145 L=L+1 GOTO 55 150 IF(IND-1) 160,155,160 155 IND=0 GOTO 50 160 IF(THR-ANRMX) 165,165,45 165 IQ=-N DO 185 I=1,N IQ=IQ+N LL=I+(I*I-I)/2 JQ=N*(I-2) DO 185 J=I,N JQ=JQ+N MM=J+(J*J-J)/2 IF(A(LL)-A(MM)) 170,185,185 170 X=A(LL) A(LL)=A(MM) A(MM)=X IF(MV-1) 175,185,175 175 DO 180 K=1,N ILR=IQ+K IMR=JQ+K X=R(ILR) R(ILR)=R(IMR) 180 R(IMR)=X 185 CONTINUE V(:)=R(:) Rx(:) = A(:) WRITE(*,*)' KET THUC EIGEN' END SUBROUTINE EIGEN !--------------------------------------------------------------------------------- SUBROUTINE TRACE(M,Nmode,K) Integer :: M, K, I, L Real :: CON WRITE(*,*)' BAT DAU C.T TRACE' R(:) = Rx(:) FM=M L=0 DO I=1,M L=L+I D(I)=R(L) ENDDO K=0 DO I=1,M !IF (D(I) < CON) Exit K=K+1 if (K>Nmode) Exit D(I)=D(I)/FM ENDDO DO I=2,K D(I)=D(I)+D(I-1) ENDDO print*,"Calculating TRACE ... Done !" END SUBROUTINE TRACE !--------------------------------------------------------------------------------- SUBROUTINE LOAD(M,K) Real :: SQ Integer :: M, K, I, J, L R(:)=Rx(:) L=0 JJ=0 DO J=1,K JJ=JJ+J SQ=SQRT(R(JJ)) DO I=1,M L=L+1 V(L)=SQ*V(L) ENDDO ENDDO print*,"Calculating LOAD ... Done !" END SUBROUTINE LOAD !----------------------------------------------------------------------------------! Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" !--------------------------------------------------------------------------------- END Program