! Chuong trinh tinh PCs, EOFs theo cach moi (2019-10-10) Use NetCDF Implicit None Real, Parameter :: Miss=-99.0 Integer, Parameter :: NMAX=100, MMAX=300 ! Cac bien Output: Real :: TriRieng(MMax),VTORieng(MMax,MMax),TaiTrong(MMax,MMax) Real :: NhanTo(NMax,MMax),TT_DAQUAY(MMax,MMax), NT_DAQUAY(NMax,MMax) Real :: XMean(MMax), Stdev(MMax), Corr(MMax,MMax) !---------------------------------------------------------------------------- Real, Allocatable :: EOFs(:,:), PCs(:,:), Eigenvalues(:), ECumulate(:) Real, Allocatable :: Mean(:), Std(:) Character (Len = 15), Allocatable :: VarName(:), Sample(:), Yname(:) Real :: XX(Nmax, Mmax) Integer :: YY(Nmax) !Real*8, Allocatable :: xdat(:,:), ydat(:,:), xtb(:), sx(:) !Real*8, Allocatable :: EOFs(:,:), PCs(:,:) !Real*8, Allocatable :: Eigenvalues(:), ECumulate(:) !--------------------------------------------------------------------------------- Integer :: Nt, Nx ! Kich thuoc thuc cua mang X Integer :: Nmode, K Character (Len=250) :: Fname, Fname1, MaskFname, pca_out, eof_out Integer :: i,j,L ! Mo ta thong tin SL dau vao: !--------------------------------------------------------------------------------- Real, allocatable :: xlon(:), xlat(:), xtime(:) Integer :: Nlon, Nlat, Ntime Logical :: Status Real, allocatable :: TC_tr(:,:,:), TC_Track(:,:) Integer :: Nst, iidx !--------------------------------------------------------------------------------- ! Reading data: Fname = "/home/tanpv/Data/TCs/WP_Annual_TC_NTC_grid_2x2.nc" Call Read_TC_track (Fname) Nx = Nst; Nmode=5 !--------------------------------------------------------------------------------- !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 (EOFs(Nt,Nmode), PCs(Nx,Nmode), ECumulate(Nx), Eigenvalues(Nx)) !--------------------------------------------------------------------------------- !Xdat(1:Nt,1:Nx) = TC_Track(1:Nt,1:Nx) ! Tinh cac EOFs, PCs !Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, & ! Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode) Call PCA_TC pca_out="fa_pc.nc"; eof_out="fa_eof.nc" Call Write_EOF_PCA (pca_out,eof_out) print*, "===========> Complete " !Deallocate(Xbar,Std,Rx,EValues,V,Fx,Ax,D,EVtor,Lamda,EOF,PCs,ECumulate) print*,"Program is normal terminating ..." CONTAINS !--------------------------------------------------------------------------------- Subroutine PCA_TC Integer :: Nth0, iyr, i, j Real :: Con CON = 1. Do iyr=1,Nt !write(Sample(iyr),'(I4)') YY(iyr) XX(iyr,1:Nx) = TC_Track(iyr,1:Nx) Enddo Nth0 = Nmode print*, Nt, Nx, Nth0 Call PTNT(XX,NMAX,MMAX,Nt,Nx,CON, Nth0, & XMean,Stdev,Corr,Nmode,TriRieng,VTORieng,TaiTrong,NhanTo,TT_DAQUAY,NT_DAQUAY) ! For other models use: EOFs(1:Nt,1:Nmode) = NT_DAQUAY(1:Nt,1:Nmode) ! NhanTo(1:Nt,1:Nmode) ! For Forecast Models: PCs (1:Nx,1:Nmode) = TT_DAQUAY(1:Nx,1:Nmode) !TaiTrong(1:Nx,1:Nmode) Eigenvalues(1:Nx) = TriRieng(1:Nx) End Subroutine PCA_TC !================================================================================= Subroutine Write_EOF_PCA (F_pca, F_eof) Character (Len=*) :: F_eof, F_pca Character (Len=255) :: var_ln, t_unit Integer :: i,j,k,iidx,it Real :: tmp Real, Allocatable :: PCA(:,:,:), xmode(:), EOF1(:,:), X1(:), Xtmp(:) Allocate (PCA(Nlon,Nlat,Nmode), xmode(Nmode), Xtmp(Nmode)) !, EOF1(Nmode,Nt)) Allocate (EOF1(Nt,Nmode)) !------------------------------------------------------------------------- ! Write to PCs file: PCA(:,:,:)=Miss; k=0 Do j=1, Nlat Do i=1, Nlon k=k+1; PCA(i,j,1:Nmode) = PCs(k,1:Nmode) Enddo Enddo Do i=1, Nmode; xmode(i)=i; Enddo ! PCA's modes var_ln="Spatial distribution of PCAs of TC tracks" t_unit="years since 1951-01-01 " Call wr_ncdf_3d (F_pca,xlon,xlat,xmode,nlon,nlat,Nmode, PCA, & "pc", trim(t_unit),trim(var_ln), & "-", "lon", "lat","mode", Miss) !------------------------------------------------------------------------- Xtmp(1:Nmode) = Eigenvalues(1:Nmode) Call wr_ncdf_1d_apd (F_pca, Xtmp, "eig", Nmode, & "Eigen values of ", "-", "mode", Miss) tmp=0.0 Do i=1,Nx !Mcount tmp=tmp+Eigenvalues(i) if (I<=Nmode) then ECumulate(i)=0 Do j=1,I ECumulate(i)=ECumulate(i)+Eigenvalues(j) Enddo endif Enddo Do i=1,Nmode ECumulate(i)=ECumulate(i)/tmp; Xtmp(i)=Eigenvalues(i)/tmp Enddo Call wr_ncdf_1d_apd (F_pca, Xtmp, "var",Nmode, & "Percent Contribution ", "-", "mode", Miss) Xtmp(1:Nmode) = ECumulate(1:Nmode) Call wr_ncdf_1d_apd (F_pca, Xtmp, "cum",Nmode, & "Cumuluted Eigenvalues", "-", "mode", Miss ) !------------------------------------------------------------------------- ! Write to EOF file: !if (Allocated(xtime)) Deallocate (xtime); Allocate (xtime(1:Nt)) !Do it=1, Nt; xtime(it)=it; Enddo ! Months !Do i=1, Nmode; EOF1(i,1:Nt)=EOFs(1:Nt,i); Enddo Do i=1, Nmode; EOF1(1:Nt,i)=EOFs(1:Nt,i); Enddo var_ln="Timeseries of eofs" Call wr_ncdf_2d (F_eof,xtime,xmode, Nt,Nmode, EOF1, "eof", & trim(var_ln), "-", "year","mode", Miss) print*,"Write_EOF_PCA: Done !!!" End Subroutine Write_EOF_PCA !--------------------------------------------------------------------------------- Subroutine Read_TC_track (Fname) Integer :: ist, ilon, ilat, it, imon, imon1, iyr Character (Len=*) :: Fname Real, Allocatable :: X0(:), X1(:,:,:), X2(:,:,:), X3(:,:,:) !Real, Allocatable :: xlon(:), xlat(:), xtime(:) Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, TC_tr, "tc") Nst=Nlon*Nlat Allocate (TC_Track(Nt,Nst)) ist=0 Do ilat=1, Nlat Do ilon=1,Nlon ist=ist+1 TC_Track(1:Nt,ist)=TC_tr(ilon,ilat,1:Nt) Enddo Enddo print*,"Nt,Nx=",Nt, Nst End Subroutine Read_TC_track !--------------------------------------------------------------------------------- Logical Function Empty_Ar1(X,N) Integer :: N,i Real :: X(N) Logical :: Empty Empty=.False. Do i=1,N if (X(i)==Miss.OR.IsNaN(X(i))) then Empty=.True.; Exit endif Enddo Empty_Ar1=Empty End Function Empty_Ar1 !----------------------------------------------------------------------------------! 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 !----------------------------------------------------------------------------------! Include "/home/tanpv/Mylibs/Factor_PCA_EOF.f90" !----------------------------------------------------------------------------------!