Program FA_ENSO_Monthly Use NetCDF Real, Parameter :: Miss=-99. Integer, Parameter :: Nidx=3 ! PDSI, SPEI12, SPI12 Character (Len=5), Parameter :: IDX_nm(Nidx)= & (/"PDSI","SPI12","SPE12"/) !---------------------------------------------------------------------------- ! Re-define: Integer :: NMAX, MMAX Real, Allocatable :: TriRieng(:),VTORieng(:,:),TaiTrong(:,:),NhanTo(:,:),TT_DAQUAY(:,:),NT_DAQUAY(:,:) Real, Allocatable :: XMean(:),Stdev(:),Corr(:,:) Real, Allocatable :: XX(:,:) Integer, Allocatable :: YY(:) !---------------------------------------------------------------------------- ! Origin Data: Real, Allocatable :: All_IDX(:,:,:,:) ! Lon,Lat,Ntime,Nidx Real, Allocatable :: Dr_IDX(:,:,:) ! Nst,Ntime,Nidx Real, Allocatable :: xlon(:), xlat(:), xtime(:), Mask(:,:) !---------------------------------------------------------------------------- Real, Allocatable :: EOFs(:,:,:), PCs(:,:,:), Eig(:,:), Var_Accum(:,:) ! ulate Real, Allocatable :: Mean(:), Std(:), Lon(:), Lat(:) Character (Len = 15), Allocatable :: VarName(:), Sample(:), Yname(:) Integer :: Nth0 ! So gia tri rieng duoc giu lai !---------------------------------------------------------------------------- Integer :: Nst, IO, yr, Ntr, iyr, ist, imon, K_Eigen Real :: Con Character (Len=4) :: Var Character (Len = 255) :: Fname, Fname1, Out_Dir, Out_Dir_Fcs Integer :: Nmon, NtrMx, NEigenMax, NTC_src Integer :: Nyear, Ntime Integer :: Nlon, Nlat, Nx1, Ny1, Nxy Character (Len=100) :: PCA_out_Fname !------------------------------------------------------------------------------------ NEigenMax = 5 ! test only Fname="/home/tanpv/Progs/Dr_SEA/Results/Dr_PDSI_SPEI_SPI.nc" Call Read_Drough_IDX (Fname) !------------------ Important: ---------------------------------------------- Mmax=Ntime; Nmax=Nst ! Dao chieu X thanh chieu T <<<<<<<<<<<<<<<<< For testing print*, Ntime, Nmax, Mmax print*, xLon print*, xLat pause Allocate (TriRieng(MMax),VTORieng(MMax,MMax),TaiTrong(MMax,MMax),NhanTo(NMax,MMax)) !,TT_DAQUAY(MMax,MMax), NT_DAQUAY(NMax,MMax)) print*,"OK1" Allocate (XMean(MMax), Stdev(MMax), Corr(MMax,MMax)) !Allocate (XMean(NMax), Stdev(NMax), Corr(MMax,MMax)) Allocate (XX(Nmax, Mmax),YY(Nmax)) !---------------------------------------------------------------------------- !Allocate (EOFs(1:Nst,1:NEigenMax,1:Nidx)) !---------------------------------------------------------------------------- !print*, GPI; pause Do iidx=1, Nidx -1 Call Calc_PCA_EOF (iidx) print*, "PCA Calculation for ",IDX_nm(iidx)," ... Done!" Enddo Fname="PCAs.nc"; Fname1="EOFs.nc" Call Write_PCA_EOF (trim(Fname),trim(Fname1)) print*, "All_Timeseries... Done!" !----------------------------------------------------------------------- CONTAINS !----------------------------------------------------------------------------------! Subroutine Calc_PCA_EOF (iidx) Integer :: ist, i, j, iidx Real :: Eigen_Acc ! Tong cac Eigen values CON = 1. !----------------- cac bien ------------ if (.not.Allocated (Sample)) Allocate (Sample(Nst)) print*,"Calc_PCA_EOF of ",IDX_nm(iidx),": ",Ntime,Nst Do ist=1, Nst ! Nyear write(Sample(ist),'(I5)') ist ! YY(iyr) !print*,"Station/Grid-point: ",ist XX(ist,1:Ntime) = Dr_IDX(ist,1:Ntime,iidx) Enddo Nth0 = NEigenMax print*, Ntime,Nst, Nth0 Call PTNT(XX,Nmax,Mmax,Nst,Ntime,CON, Nth0, & XMean,Stdev,Corr,K_Eigen,TriRieng,VTORieng,TaiTrong,NhanTo) !,TT_DAQUAY,NT_DAQUAY) print*,"Done PCA_EOF calling.",Nst, Ntime print*,"EOF shape: ",shape(NhanTo) print*,"PCA shape: ",shape(TaiTrong) print*,"Nth0=", Nth0,NEigenMax !pause NEigenMax = Nth0 ! For other models use: if (.not.Allocated(EOFs)) Allocate (EOFs(1:Nst,1:NEigenMax,1:Nidx)) if (.not.Allocated(PCs)) Allocate (PCs(1:Ntime,1:NEigenMax,1:Nidx)) if (.not.Allocated(Eig)) Allocate (Eig(1:NEigenMax,1:Nidx)) if (.not.Allocated(Var_Accum)) Allocate (Var_Accum(1:NEigenMax,1:Nidx)) EOFs(1:Nst,1:NEigenMax,iidx) = NhanTo(1:Nst,1:NEigenMax) PCs (1:Ntime,1:NEigenMax,iidx) = TaiTrong(1:Ntime,1:NEigenMax) Eig(1:NEigenMax,iidx) = TriRieng(1:NEigenMax) Eigen_Acc=0 Do i=1, Mmax Eigen_Acc=Eigen_Acc+TriRieng(i) Enddo Do i=1, NEigenMax Var_Accum(i,iidx)=0 Do j=i,NEigenMax Var_Accum(i,iidx)=Var_Accum(i,iidx)+TriRieng(j) Enddo Var_Accum(i,iidx)=Var_Accum(i,iidx)/Eigen_Acc*100.0 Enddo !Mean(1:Ntime) = XMean(1:Ntime) !Std(1:Ntime) = Stdev(1:Ntime) !print*,"<<<<<<<<<<<<<<" !----------------------------------------------------------------------- ! Testing: ! print*,"FT x F = Lamda" ! print "(10F8.2)", matmul (transpose(PCs),PCs); pause ! print*,"AT x A = I" ! print "(10F12.2)", matmul (transpose(EOFs),EOFs); pause !! Enddo ! k=1,Nmon End Subroutine Calc_PCA_EOF !----------------------------------------------------------------------- Subroutine Write_PCA_EOF (Fnm_PCA, Fnm_EOF) Character (Len=*) :: Fnm_PCA, Fnm_EOF Character (Len=100) :: Var_nm, Var_ln, t_unit Real, Allocatable :: X(:,:,:), Y(:,:), Y1(:), Y2(:), Ymean(:,:), Ystd(:,:), xmod(:) Integer :: ilon,ilat,ik,it, ist, iidx Allocate (X(1:Nlon,1:Nlat,1:NEigenMax),xmod(NEigenMax))! ,Y(1:Nlon,1:Nlat)) Do ik=1, NEigenMax; Xmod(ik)=ik; Enddo t_unit="months since 1977-01-01 " Do iidx=1, Nidx ist=0; X(:,:,:)=Miss Do ilat=1, Nlat Do ilon=1, Nlon if (Mask(ilon,ilat) >0) then ist=ist+1; X(ilon,ilat,1:NEigenMax)=EOFs(ist,1:NEigenMax,iidx) endif Enddo Enddo Var_nm=IDX_nm(iidx) Var_ln="PCA of "//IDX_nm(iidx)//" drought index" if (iidx==1) then Call wr_ncdf_3d (Fnm_PCA,xlon,xlat,Xmod,nlon,nlat,NEigenMax,X, trim(Var_nm), trim(t_unit), & trim(Var_ln), "-", "lon", "lat","mod", Miss) else Call wr_ncdf_3d_apd (Fnm_PCA, X, trim(var_nm), nlon,nlat,NEigenMax, trim(Var_ln), & "-", "lon", "lat","mod", Miss) endif !Var_nm="Eig_"//trim(IDX_nm(iidx)) !Var_ln="Eigen values of "//IDX_nm(iidx)//" drought index" !Call wr_ncdf_1d_apd (Fnm_PCA, Eig, trim(Var_nm), NEigenMax,trim(Var_ln), "-", "mod", Miss) Enddo Deallocate (X) !Allocate (Ymean(1:Nlon,1:Nlat), Ystd(1:Nlon,1:Nlat)) !it=0 !Do ilat=1, Nlat ! Do ilon=1, Nlon ! it=it+1; Ymean(ilon,ilat)=Mean(it); Ystd(ilon,ilat)=Std(it) ! Enddo !Enddo !Call wr_ncdf_2d_apd (Fnm_PCA, Ymean, "mean",nlon,nlat,"GPI mean", "-", "lon", "lat", Miss) !Call wr_ncdf_2d_apd (Fnm_PCA, Ystd, "std" ,nlon,nlat,"GPI_Std", "-", "lon", "lat", Miss) !Call wr_ncdf_1d_apd (Fnm_PCA, Eig, "eigen",NEigenMax,"Eigen values", "-", "mod", Miss) !Deallocate (Ymean,Ystd) Allocate (Y(1:Ntime,1:NEigenMax), Y1(1:NEigenMax), Y2(1:NEigenMax)) Do iidx=1, Nidx Do ik=1, NEigenMax Y(1:Ntime,ik) = PCs (1:Ntime,ik,iidx) Y1(ik) = Eig(ik,iidx) Y2(ik) = Var_Accum(ik,iidx) Enddo Var_nm=IDX_nm(iidx) Var_ln="Timeseries of EOFs of "//IDX_nm(iidx)//" drought index" if (iidx==1) then Call wr_ncdf_2d (Fnm_EOF, Xtime,Xmod, Ntime,NEigenMax, Y, trim(var_nm), trim(Var_ln), "-", "mon", "mod", Miss) else Call wr_ncdf_2d_apd (Fnm_EOF, Y, trim(var_nm), Ntime,NEigenMax, trim(Var_ln), "-", "mon", "mod", Miss) endif Var_nm="Eig_"//trim(IDX_nm(iidx)) Var_ln="Eigen values of "//IDX_nm(iidx)//" drought index" Call wr_ncdf_1d_apd (Fnm_EOF, Y1, trim(var_nm), NEigenMax, trim(Var_ln), "-", "mod", Miss) Var_nm="Acu_Eig_"//trim(IDX_nm(iidx)) Var_ln="Accumulated Eigen values of "//IDX_nm(iidx)//" drought index" Call wr_ncdf_1d_apd (Fnm_EOF, Y2, trim(var_nm), NEigenMax, trim(Var_ln), "-", "mod", Miss) Enddo End Subroutine Write_PCA_EOF !------------------------------------------------------------------------------------ Subroutine Read_Drough_IDX (Fname) Character (Len=8), Parameter :: Vname(Nidx)=(/"ori_pdsi","SPI_12","SPEI_12"/) Integer :: iidx, ist, ilon, ilat, it, Nt, imon, imon1, iyr Integer :: Nyear, Nmon, nlat1, nlat0, dlat Character (Len=*) :: Fname Real :: tmp, xmean(12), ic(12) Real, Allocatable :: X0(:), X(:,:), X1(:,:,:), X2(:,:,:), X3(:,:,:), XDr(:,:,:) !PDSI(:,:,:), SPI12(:,:,:), SPEI12(:,:,:), Mask(:,:) !xlon(:), xlat(:), xtime(:), Nlon, Nlat, Ntime Nyear=30 ! 1976-2005 Nmon =12 dlat = 16 ! So diem vi do se bo (o vung Australia) Do iidx=1, Nidx if (Allocated(xlon)) Deallocate (xlon) if (Allocated(xlat)) Deallocate (xlat) if (Allocated(xtime)) Deallocate (xtime) Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, XDr, trim(Vname(iidx))) if (.not.Allocated (All_IDX)) Allocate (All_IDX(nlon,nlat,Nt-Nmon,Nidx)) All_IDX(:,:,1:Nt-Nmon,iidx) = XDr(:,:,Nmon+1:Nt); Deallocate (XDr) Enddo Deallocate (xlon, xlat) Call rd_ncdf_2d (trim(Fname), xlon,xlat, nlon,nlat, X, "land") ! Mask,"land") nlat1 = nlat - dlat Ntime = Nt - Nmon Allocate (Mask(nlon,nlat1), lat(nlat1)) Mask(1:nlon,1:nlat1) = X(1:nlon,dlat+1:nlat) ! Cho mien da bo vung Australia lat(1:nlat1) = xlat(dlat+1:nlat) Deallocate (xlat); Allocate (xlat(Nlat1)) xlat(1:Nlat1) = lat(1:nlat1) print*,"Read_Drough_IDX: ... done !" ist=0 Do ilat=1, Nlat1 Do ilon=1, Nlon if (Mask(ilon,ilat)>0) ist=ist+1 Enddo Enddo Nst=ist Allocate (Dr_IDX(Nst,Ntime,Nidx)) Do iidx=1, Nidx ist=0 Do ilat=1, Nlat1 Do ilon=1, Nlon if (Mask(ilon,ilat)>0) then ist=ist+1 Dr_IDX(ist,1:Ntime,iidx)=All_IDX(ilon,ilat+dlat,1:Ntime,iidx) endif Enddo Enddo Enddo nlat=nlat1 print*,"==========" End Subroutine Read_Drough_IDX !----------------------------------------------------------------------- Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" END ! Main Prog !----------------------------------------------------------------------- !Include "Old/Factor_PCA_EOF.f90" Include "/home/tanpv/Mylibs/Sub_PCA_EOF_New.f90" !Include "/home/tanpv/Mylibs/Factor_PCA_EOF.f90"