Program FA_ENSO_Monthly ! Chuong trinh nay thuc hien viec so sanh ket qua PCA theo hai phuong an: ! PA1: Tinh theo bo chuong trinh cu (/home/tanpv/Mylibs/Factor_PCA_EOF.f90) ! PA2: Tinh theo bo chuong trinh moi (/home/tanpv/Mylibs/Sub_PCA_2019.f90) ! Viec tinh toan duoc thuc hien tren cung mot bo SL ClimIDX ! Ket qua PA1 ghi ra file FA_ClimIDX.txt ! Ket qua PA2 ghi ra file EOF_New_ClimIDX.txt ! Viec so sanh noi dung hai file cho thay: ! PA1 giong voi PA2 khi PA2 chon cac tham so: ! switch=1; Standard=1; PC_Stard=0 ! hoac switch=0; Standard=1; PC_Stard=0 Use NetCDF Integer, Parameter :: NMAX=1000, MMAX=20 Real, Parameter :: Miss=-99. ! Cac bien Output: Real :: TriRieng(MMax),VTORieng(MMax,MMax),TaiTrong(MMax,MMax),NhanTo(NMax,MMax),TT_DAQUAY(MMax,MMax), NT_DAQUAY(NMax,MMax) Real :: XMean(MMax), Stdev(MMax), Corr(MMax,MMax) Real :: XX(Nmax, Mmax) Integer :: YY(Nmax) !---------------------------------------------------------------------------- ! Origin Data: Real, Allocatable :: ClimIDX(:,:,:) Real*8, Allocatable :: xdat(:,:), ydat(:,:), xtb(:), sx(:) Real*8, Allocatable :: EOFss(:,:), PCss(:,:) Real*8, Allocatable :: Eigenvalues(:), ECumulate(:) Integer :: Nt, Nx ! Kich thuoc thuc cua mang X Integer :: Mcount, Ncount, switch, Standard, PC_Stard, Nmode, K Character (Len=250) :: Fname, Fname1, MaskFname, pca_out, eof_out Real*8 :: Factor !---------------------------------------------------------------------------- Real, Allocatable :: EOFs(:,:), PCs(:,:), Eig(:) Real, Allocatable :: Mean(:), Std(:) Character (Len = 15), Allocatable :: VarName(:), Sample(:), Yname(:) Integer :: Nth0 ! So gia tri rieng duoc giu lai !---------------------------------------------------------------------------- Integer :: NTram, IO, yr, Ntr, iyr, ist, imon, K_Eigen Real :: Con Character (Len=4) :: Var Character (Len = 255) :: Out_Dir, Out_Dir_Fcs Integer :: NYr1, NYr2, Nmon, NtrMx, NEigenMax, NTC_src Integer :: Nyear, Ntime Integer :: Kmon, Lag Integer :: mon_y Integer :: Nyrs, yrs_y Character (Len=100) :: ClimFname, TCFname, PCA_out_Fname Namelist /Data_info_Training/ NYr1, NYr2, Nmon, NtrMx, NEigenMax Namelist /ClimIDX_Name/ VarName Namelist /PCA_inp_Out_Fname/ ClimFname !---------------------------------------------------------------------------- Open (1, file="NameList_FA.txt", Status="Old") Read (1, Data_info_Training); print*, NYr1, NYr2, Nmon, NtrMx, NEigenMax Nyear = (Nyr2-Nyr1+1); Ntime = Nyear*Nmon; print*, Nyear, Ntime Allocate (VarName(NtrMx), Sample(Nmax), Yname(NTC_src)) Allocate (ClimIDX(1:Nyear,1:Nmon,1:NtrMx)) Allocate (EOFs(1:Ntime,1:NEigenMax), PCs(1:NtrMx,1:NEigenMax), Eig(1:NEigenMax)) Allocate (Mean(1:NtrMx), Std(1:NtrMx)) Read (1, ClimIDX_Name) Read (1, PCA_inp_Out_Fname) Close (1) Ntr = NtrMx; !Out_Dir="Out/" ClimFname, TCFname, PCA_out_Fname Call Read_Clim_IDX(trim(ClimFname)) ; print*, "Call Read_Clim_IDX ... Done!" Call PCA_for_Monthly_Data !(Fname) Call Write_PCA_to_TXT ("FA_ClimIDX.txt",NEigenMax,Ntime,Ntr) print*, "PCA Monthly... Done!" !----------------------------------------------------------------------- Nt = Ntime; Nx = Ntr switch=0; Standard=1; PC_Stard=0; Nmode=NEigenMax !--------------------------------------------------------------------------------- if (switch==1) then Mcount=Nt; Ncount=Nx else if (switch==0) then Mcount=Nx; Ncount=Nt endif Allocate (xdat(Nt,Nx), ydat(Nt,Nx)) Allocate (Xtb(Nx), Sx(Nx)) ! TB, DLC, ma tran Covar hoac Corre Allocate (Eigenvalues(Mcount)) ! Tri rieng va Vecto rieng (1D) Allocate (EOFss(Nt,Nmode), PCss(Nx,Nmode), ECumulate(Nmode)) !--------------------------------------------------------------------------------- xdat(1:Nt,1:Nx)=XX(1:Nt,1:Nx) ! Tinh cac EOFs, PCs Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, & Xtb,Sx,Eigenvalues,EOFss,PCss,Nmode) Deallocate(EOFs,PCs,Eig) Allocate (EOFs(1:Ntime,1:NEigenMax), PCs(1:NtrMx,1:NEigenMax), Eig(1:NEigenMax)) EOFs(1:Ntime,1:Nmode)=EOFss(1:Ntime,1:Nmode) PCs(1:Nx,1:Nmode)=PCss(1:Nx,1:Nmode) Eig(1:Nmode)=Eigenvalues(1:Nmode) Call Write_PCA_to_TXT ("EOF_New_ClimIDX.txt",NEigenMax,Ntime,Ntr) CONTAINS !----------------------------------------------------------------------- Subroutine Read_Clim_IDX (Fname) ! Read data of Climate Indices from NetCDF and store in 3D array of ClimIDX(1:Nyear,1:Nmon,1:Ntr) ! where Ntr=14 is number of Climate indices Integer :: iidx, iyr, imon, year, IO, it, it0, it1 Character (Len=*) :: Fname Real, Allocatable :: X(:,:),Xyear(:),Xmon(:) Integer :: Nyr, Nmn Do iidx=1, Ntr if (Allocated(X)) Deallocate (X) if (Allocated(Xyear)) Deallocate (Xyear) if (Allocated(Xmon)) Deallocate (Xmon) Call rd_ncdf_2d(trim(Fname),Xyear,Xmon,Nyr,Nmn, X, trim(VarName(iidx))) Do iyr=NYr1, NYr2 Do it=1, Nyr if (Int(Xyear(it))==iyr) then it0=iyr-NYr1+1 ClimIDX(it0,:,iidx)=X(it,:) YY(it0) = iyr endif Enddo Enddo ! Create XX array for FA: it=0 Do iyr=1,Nyear Do imon=1,Nmon it=it+1; XX(it,iidx)=ClimIDX(iyr,imon,iidx) Enddo Enddo Enddo End Subroutine Read_Clim_IDX !----------------------------------------------------------------------- Subroutine PCA_for_Monthly_Data !(Fname) ! Calculate the PCA and EOF and store in the 2D arrays of TaiTrong, NhanTo, respectively ! for each month of the year ! EOFs(1:Nyear,1:NEigenMax,1:Nmon) !Character (Len=*) :: Fname !Character (Len=255) :: Fname1, Fmt Integer :: iyr, imon, i, j CON = 1. !----------------- cac bien ------------ i=0 Do iyr=NYr1, NYr2 Do imon=1, Nmon i=i+1 write(Sample(i),'(I4.4,I2.2)') iyr, imon Enddo Enddo Nth0 = NEigenMax print*, Ntime, Ntr, Nth0 Call PTNT(XX,NMAX,MMAX,Ntime,Ntr,CON, Nth0, & XMean,Stdev,Corr,K_Eigen,TriRieng,VTORieng,TaiTrong,NhanTo,TT_DAQUAY,NT_DAQUAY) ! For other models use: EOFs(1:Ntime,1:NEigenMax) = NhanTo(1:Ntime,1:NEigenMax) ! For Forecast Models: PCs (1:NtrMx,1:NEigenMax) = TaiTrong(1:Ntr,1:NEigenMax) Eig(1:NEigenMax) = TriRieng(1:NEigenMax) Mean(1:NtrMx) = XMean(1:Ntr) Std(1:NtrMx) = Stdev(1:Ntr) End Subroutine PCA_for_Monthly_Data !----------------------------------------------------------------------- Subroutine Write_PCA_to_TXT (Fname,K_Eigen,N_Sample,M_Var) Integer :: K_Eigen,N_Sample,M_Var Character (Len=*) :: Fname Open(4,file=trim(Fname), status="Unknown") write(4,"(I3,A,16F8.3)") K_Eigen, " tri rieng: ",(Eig(i),i=1,K_Eigen) ! (TriRieng(i),i=1,K_Eigen) !write(4,*) "Cac vector rieng:" !Do j=1,K_Eigen ! write(4,"(I3,A,16F8.3)") j," : ", (VTORieng(i,j),i=1,M_Var) !Enddo !write(4,*) "Ma tran Tai trong nhan to (PCA):" write(4,*) "Ma tran cac PCs" Do i=1, M_Var write(4,"(A11,16F8.4)") VarName(i), (PCs(i,j),j=1,K_Eigen) !write(4,"(A11,16F8.4)") VarName(i), (TaiTrong(i,j),j=1,K_Eigen) Enddo !write(4,*) "Ma tran Tai trong nhan to da quay (PCA):" !Do i=1, M_Var ! write(4,"(A11,16F8.4)") VarName(i), (TT_DAQUAY(i,j),j=1,K_Eigen) !Enddo !write(4,*) "Ma tran nhan to (EOF):" write(4,*) "Ma tran cac EOFs" Do i=1, N_Sample write(4,"(A11,16F8.4)") Sample(i), (EOFs(i,j),j=1,K_Eigen) !(NhanTo(i,j),j=1,K_Eigen) Enddo !write(4,*) "Ma tran nhan to da quay (EOF):" !Do i=1, N_Sample ! write(4,"(A11,16F8.4)") Sample(i), (NT_DAQUAY(i,j),j=1,K_Eigen) !Enddo Close(4) End Subroutine Write_PCA_to_TXT !----------------------------------------------------------------------- Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Sub_PCA_2019.f90" END ! Main Prog !----------------------------------------------------------------------- !Include "Old/Factor_PCA_EOF.f90" Include "/home/tanpv/Mylibs/Factor_PCA_EOF.f90"