! Chuong trinh tinh PCs, EOFs theo cach moi (2019-10-10) Use NetCDF Implicit None Real, Parameter :: Miss=-99.0, Lat_min=-11.0 ! Cat phan Autralia Integer, Parameter :: Nmon=12, Nidx=16 !24 Integer, Parameter :: Nlen=3, Ndr_idx=2 ! (3,6,12 month); spi and spei Integer, Parameter :: Len_SPI(Nlen)=(/3,6,12/) Character (len=4), parameter :: idx_nm(Ndr_idx)=(/"spi","spei"/) Character (len=2), parameter :: Timescale(Nlen)=(/"03","06","12"/) Real*8, Allocatable :: xdat(:,:), ydat(:,:), xtb(:), sx(:) Real*8, Allocatable :: EOFs(:,:), PCs(:,:) Real*8, Allocatable :: Eigenvalues(:), ECumulate(:) !--------------------------------------------------------------------------------- Real, Allocatable :: PC_GrADS(:,:,:) ! (Nst,Nmode,Nlen) Save all PCs for GrADS Real, Allocatable :: PC_all(:,:,:), Eof_all(:,:,:), Eig_all(:,:) ! PCs,Eigenvalues,EOFs !--------------------------------------------------------------------------------- Integer :: Yr1, Yr2, Nyear ! Thoi ky Integer :: Nt, Nx ! Kich thuoc thuc cua mang X Integer :: Mcount, Ncount, switch, Standard, PC_Stard, Nmode, K Integer :: switch1, standard1, PC_Stard1, Nmode1 ! Clim Integer :: switch2, standard2, PC_Stard2, Nmode2 ! sst Character (Len=250) :: Fname, Fname1, MaskFname, pca_out, eof_out Character (Len=250) :: Fnclim, pca_clim, eof_clim Character (Len=250) :: Fn_dr, pca_dr, eof_dr Real :: Slon1, Slon2, Slat1, Slat2, Miss_sst ! sst domain Integer :: Yr_sta Real*8 :: Factor Namelist /Model_Info/ Yr1, Yr2, switch, standard, PC_Stard, Nmode Namelist /Clim_info/ Fnclim, pca_clim, eof_clim, & Fn_dr, pca_dr, eof_dr !switch1, standard1, PC_Stard1, Nmode1 Integer :: i,j,L ! Mo ta thong tin SL dau vao: !--------------------------------------------------------------------------------- Real, allocatable :: xlon(:), xlat(:), Reg(:), xtime(:) Character (Len=20), Allocatable :: stname(:) Integer :: Ntime, ist Logical :: Status Real, allocatable :: SPI(:,:,:), SPEI(:,:,:) Real, allocatable :: ClimIDX(:,:) Real, allocatable :: Dr_IDX(:,:,:,:) ! (Nst,Ntime,Nlen,NDr_idx) Real, allocatable :: Mask(:,:), Mask_all(:,:) !--------------------------------------------------------------------------------- ! SST: Real, allocatable :: lon_sst(:), lat_sst(:) Real, allocatable :: sst(:,:), Mask_sst(:,:) ! (Nt,Nx) Integer :: Nlon_sst, Nlat_sst, SeaPoint, Mov_Ave_Step !--------------------------------------------------------------------------------- ! Composite maps of SST: Real, allocatable :: sst_dry(:,:), sst_wet(:,:) ! (Nlon_sst, Nlat_sst): SST with Areal mean of SPEI<-0.5 and >-0.5 Real, allocatable :: sst_3d(:,:,:) ! (Nlon_sst, Nlat_sst,Ntime) Real, allocatable :: Areal_mean_Dr(:) ! (Ntime) !--------------------------------------------------------------------------------- ! Rainfall: Real, allocatable :: Prec(:,:), lon_prec(:), lat_prec(:), Mask_prec(:,:) Integer :: Ntime_prec, Nlon_prec, Nlat_prec, PrecPoint Integer :: Nst, iidx, ilen, it !Integer, Parameter :: NDr_idx=2 !Character (Len=5), Parameter :: IDX_nm(NDr_idx)=(/spi","spei"/) ! (/spi","spei"/) Character (Len=255) :: Filename !--------------------------------------------------------------------------------- Character (Len=15), Parameter :: Clim_nm(Nidx)=(/ & ! "nino3", "nino34", "nino4", "oni", "censo", & ! "tni", & ! 6 ! "rindo_slpa","soi", "reqsoi", "repac_slpa", & !"whwp", "ammsst", "pdo", "NTA_ersst", "dmi", & ! 5 Thay CAR_ersst bang pdo, thay tsa bang dmi ! "whwp", "ammsst", "NTA_ersst", "tna", "amonus", & ! "pdo", "dmi"/) ! "pna", "censo", "soi", "reqsoi", "np", "repac_slpa", & ! 6 "nino12", "nino3", "nino34", "nino4", "oni", "tni", & ! 6 "whwp", "pdo", "dmi", "tna", "amonus", & "pna", "censo", "soi", "reqsoi", "repac_slpa"/)!, & ! 6 ! "wp"/) ! "nino12", "nino3", "nino34", "nino4", "oni", "tni", & ! 6 ! "whwp", "ammsst", "pdo", "NTA_ersst", "dmi", & ! "tna", "amonus", & ! "pna", "censo", "soi", "reqsoi", "np", "repac_slpa", & ! 6 ! "rindo_slpa", "qbo30", "wp", "ao", "nao"/) !--------------------------------------------------------------------------------- Open (1,file="Namelist_PCA_Dr_VN_2025.txt",Status="old") Read (1, Model_Info) Read (1, Clim_info) !Read (1, SST_info) Close (1) Mov_Ave_Step = 13 ! Moving average !--------------------------------------------------------------------------------- Nyear = Yr2-Yr1 ! Bo nam dau tien Ntime = Nyear*Nmon !--------------------------------------------------------------------------------- ! PCA for Clim: !--------------------------------------------------------------------------------- Call Read_ClimIDX (Fnclim) Nt = Ntime; Nx = Nidx !switch=switch1; standard=standard1; PC_Stard=PC_Stard1; Nmode=Nmode1 PC_Stard=1 ! Test !Nmode=6 if (switch==1) then Mcount=Nt; Ncount=Nx else if (switch==0) then Mcount=Nx; Ncount=Nt endif Factor = Sqrt(Real(Nt)) 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(Nmode)) Xdat(1:Nt,1:Nx) = ClimIDX(1:Nt,1:Nx) Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, & Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode) pca_clim=trim("/home/tanpv/Progs/Dr_4Books/Res/")//trim(pca_clim) eof_clim=trim("/home/tanpv/Progs/Dr_4Books/Res/")//trim(eof_clim) Call Write_EOF_PCA_Clim (pca_clim, eof_clim) !Call Print_test; pause Deallocate (xdat,ydat,Xtb,Sx,Eigenvalues,EOFs,PCs,ECumulate) !print*,"Program is normal terminating ..." print*, "===========> Complete for ClimIDX" !--------------------------------------------------------------------------------- ! PCA for drought: !--------------------------------------------------------------------------------- Call Read_Dr_IDX (Fn_dr) Nt = Ntime; Nx = Nst !Nmode=6 if (switch==1) then Mcount=Nt; Ncount=Nx else if (switch==0) then Mcount=Nx; Ncount=Nt endif Factor = Sqrt(Real(Nt)) Do iidx=1,Ndr_idx if (Allocated (PC_GrADS)) Deallocate (PC_GrADS) if (Allocated (PC_all)) Deallocate (PC_all) if (Allocated (Eof_all)) Deallocate (Eof_all) if (Allocated (Eig_all)) Deallocate (Eig_all) Allocate (PC_GrADS(Nx,Nmode,Nlen)) Allocate (PC_all(Nx,Nmode,Nlen), Eof_all(Nt,Nmode,Nlen)) Allocate (Eig_all(Mcount,Nlen)) Do ilen=1,Nlen 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(Nmode)) Do ist=1,Nx Do it=1,Nt Xdat(it,ist) = Dr_IDX(ist,it,ilen,iidx) Enddo Enddo Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, & Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode) PC_GrADS(1:Nx,1:Nmode,ilen) = PCs (1:Nx,1:Nmode) PC_all (1:Nx,1:Nmode,ilen) = PCs (1:Nx,1:Nmode) Eof_all (1:Nt,1:Nmode,ilen) = EOFs(1:Nt,1:Nmode) Eig_all (1:Mcount,ilen) = Eigenvalues(1:Mcount) !pca_clim=trim("/home/tanpv/Progs/Dr_4Books/Res/")//trim(pca_dr)//"_"//Timescale(ilen)//".nc" !eof_clim=trim("/home/tanpv/Progs/Dr_4Books/Res/")//trim(eof_dr)//"_"//Timescale(ilen)//".nc" !Call Write_EOF_PCA_Clim (pca_clim, eof_clim) print*,idx_nm(iidx),"--",Timescale(ilen) !Call Print_test; pause Fname = trim("/home/tanpv/Progs/Dr_4Books/Res/")//trim(pca_dr)//"_"// & trim(idx_nm(iidx))//"_"//Timescale(ilen)//"_for_Clustering.nc" Call Write_PCA_for_Clustering (Fname) Deallocate (xdat,ydat,Xtb,Sx,Eigenvalues,EOFs,PCs,ECumulate) Enddo pca_clim=trim("/home/tanpv/Progs/Dr_4Books/Res/")//trim(pca_dr)//"_"//trim(idx_nm(iidx))//".nc" eof_clim=trim("/home/tanpv/Progs/Dr_4Books/Res/")//trim(eof_dr)//"_"//trim(idx_nm(iidx))//".nc" Call Write_EOF_PCA_Dr (pca_clim, eof_clim) Fname = trim("/home/tanpv/Progs/Dr_4Books/Res/")//trim(pca_dr)//"_"//trim(idx_nm(iidx)) Call Write_PCs_to_GrADS (Fname) Enddo !print*,"Program is normal terminating ..." print*, "===========> Complete for ClimIDX" !--------------------------------------------------------------------------------- CONTAINS !================================================================================= Subroutine Write_PCs_to_GrADS (Fname) Character (Len=*) :: Fname Character (Len=4),Allocatable :: vname(:) ! (/"mod1","mod2","mod3","mod4","mod5","mod6"/) Allocate (vname(Nmode)) vname(:)=(/"mod1","mod2","mod3","mod4","mod5","mod6"/) Call Write_GrADS_Stn_3D (vname,PC_GrADS,Nst,Nmode,Nlen,xlon,xlat,Stname,Fname) End Subroutine Write_PCs_to_GrADS !------------------------------------------------------------------------- Subroutine Write_PCA_for_Clustering (Fname) Character (Len=*) :: Fname Real, Allocatable :: PCA(:,:),xmode(:),xpoint(:) ! (Nx,Nmode) Integer :: i,j,k,iidx,it if (Allocated (PCA)) Deallocate (PCA) if (Allocated (xmode)) Deallocate (xmode) if (Allocated (xpoint)) Deallocate (xpoint) Allocate (PCA(Nx,Nmode),xmode(Nmode),xpoint(Nx)) Do i=1, Nmode; xmode(i)=i; Enddo Do i=1,Nx; xpoint(i)=i; Enddo Do i=1,Nmode; PCA(1:Nx,i)=PCs(1:Nx,i); Enddo !if (iidx==1.OR.iidx<1.OR.iidx>3) then Call wr_ncdf_2d (Fname,xpoint,xmode,Nx,Nmode,PCA,"pca", & "PC values of Drought index", "-", "point","mod", Miss) !else ! Call wr_ncdf_2d_apd (Fname,PCA,"???",Nx,Nmode, & ! trim(var_ln), "-", "mon", "mode", Miss) !endif End Subroutine Write_PCA_for_Clustering !------------------------------------------------------------------------- Subroutine Write_EOF_PCA_Dr (F_pca, F_eof) ! Allocate (PC_all(Nx,Nmode,Nlen), Eof_all(Nt,Nmode,Nlen)) ! Allocate (Eig_all(Mcount,Nlen)) Character (Len=*) :: F_eof, F_pca Character (Len=255) :: var_ln,var_unit,t_unit,lon_nm,lat_nm,t_nm Integer :: i,j,k,iidx,it Real :: tmp,sumx Real, Allocatable :: PCA(:,:,:), xmode(:), EOF1(:,:,:), Xstn(:), Xtmp(:) Real, Allocatable :: Percent(:,:), Accum(:,:), xtime(:),Eigen(:,:) Real :: xlen(Nlen) Allocate (PCA(Nmode,Nlen,Nx),EOF1(Nmode,Nlen,Nt)) Allocate (xmode(Nmode), Xtmp(Nmode),Xstn(Nx),xtime(Nt)) Allocate (Percent(Nmode,Nlen),Accum(Nmode,Nlen),Eigen(Nmode,Nlen)) !------------------------------------------------------------------------- Do i=1,Nmode Do j=1,Nlen PCA (i,j,1:Nx) = PC_all (1:Nx,i,j) EOF1(i,j,1:Nt) = Eof_all(1:Nt,i,j) ENddo Enddo Do i=1,Nmode; xmode(i)=i; Enddo ! mode xlen(:) = Len_SPI(:) ! Timescale Do i=1,Nx; Xstn(i)=i; Enddo ! Station Do i=1, Nt; xtime(i)=i; Enddo ! Time ! Xtmp(1:Nmode) = Eigenvalues(1:Nmode) ! Eigen values PCs,Eigenvalues,EOFs ! Accum (%) of Eigen values values: Do j=1,Nlen Sumx=Sum(Eig_all(1:Mcount,j)); !Accum(:,:)=0 Xtmp(1:Nmode) = Eig_all(1:Nmode,j) Do i=1,Nmode Percent(i,j) = Xtmp(i)/Sumx*100 Accum (i,j) = Sum(Percent(1:i,j)) !tmp=0 !Do k=1,i ! tmp=tmp+Percent(i,j) !Enddo !Accum (i,j) = tmp Eigen (i,j) = Eig_all(i,j) Enddo Enddo var_unit="-" var_ln="PCA of ClimIDXs" t_unit="months since 1981-01-01"; lon_nm="mode"; lat_nm="timescale"; t_nm="stnid" Call wr_ncdf_3d (F_pca,xmode,xlen,Xstn,Nmode,Nlen,Nx,PCA,"pca", & t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) !var_ln="Timeseries of eofs" !Call wr_ncdf_3d_apd (F_pca,EOF1,"eof",Nmode,Nlen,Nt, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) var_ln="Percent Contribution of PCs" Call wr_ncdf_2d_apd (F_pca,Percent,"per",Nmode,Nlen, var_ln, var_unit, lon_nm, lat_nm, Miss) var_ln="Accumulated Contribution of PCs" Call wr_ncdf_2d_apd (F_pca,Accum,"cum",Nmode,Nlen, var_ln, var_unit, lon_nm, lat_nm, Miss) var_ln="Eigen values of PCs" Call wr_ncdf_2d_apd (F_pca,Eigen,"eig",Nmode,Nlen, var_ln, var_unit, lon_nm, lat_nm, Miss) var_ln="Timeseries of eofs"; t_nm="month" Call wr_ncdf_3d(F_eof,xmode,xlen,xtime,Nmode,Nlen,Nt,EOF1,"eof", & t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! Call wr_ncdf_2d (F_eof,xtime,xmode, Nt,Nmode, EOF1,"eof", & ! trim(var_ln), "-", "mon","mod", Miss) End Subroutine Write_EOF_PCA_Dr !------------------------------------------------------------------------- Subroutine Write_EOF_PCA_Clim (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,sumx Real, Allocatable :: PCA(:,:), xmode(:), EOF1(:,:), X1(:), Xtmp(:) Real, Allocatable :: Percent(:), Accum(:), xtime(:) Allocate (PCA(Nmode,Nx),xmode(Nmode), Xtmp(Nmode),EOF1(Nt,Nmode)) Allocate (Percent(Nmode),Accum(Nmode),X1(Nx),xtime(Nt)) !------------------------------------------------------------------------- Do i=1,Nmode; PCA(i,1:Nx)=PCs(1:Nx,i); Enddo Do i=1, Nmode; xmode(i)=i; Enddo Do i=1,Nx; X1(i)=i; Enddo ! Clim-idx Xtmp(1:Nmode) = Eigenvalues(1:Nmode) ! Eigen values ! Accum (%) of Eigen values values: Sumx=Sum(Eigenvalues(1:Mcount)); Accum(:)=0 Do i=1,Nmode Percent(i)=Xtmp(i)/Sumx*100 Accum(i)=Sum(Percent(1:i)) Enddo var_ln="PCA of ClimIDXs" Call wr_ncdf_2d (F_pca,xmode,X1,Nmode,Nx,PCA,"pca", & trim(var_ln), "-", "mod","idx", Miss) Call wr_ncdf_1d_apd (F_pca, Percent, & "per",Nmode, & "Percent Contribution of clim idx", & "-", "mod", Miss ) Call wr_ncdf_1d_apd (F_pca, Accum, & "cum",Nmode, & "Accumulated Contribution of clim idx", & "-", "mod", Miss ) Call wr_ncdf_1d_apd (F_pca, Xtmp, & "eig",Nmode, & "Eigen values of clim idx", & "-", "mod", Miss ) !------------------------------------------------------------------------- Do i=1, Nt; xtime(i)=i; EOF1(i,1:Nmode)=EOFs(i,1:Nmode); Enddo ! Months var_ln="Timeseries of eofs of ClimIDX" Call wr_ncdf_2d (F_eof,xtime,xmode, Nt,Nmode, EOF1,"eof", & trim(var_ln), "-", "mon","mod", Miss) End Subroutine Write_EOF_PCA_Clim !--------------------------------------------------------------------------------- Subroutine Read_ClimIDX (Fname) Character (Len=*) :: Fname Real, Allocatable :: Year(:),Month(:), tmp(:,:) Integer :: Yr11, Nyr, Nmn, iidx, iyr, i,j, k Yr11=Yr1+1 ! Bo nam dau tien Ntime=(Yr2-Yr11+1)*Nmon Allocate (ClimIDX(Ntime,Nidx)) !print*,Nidx,Ntime,Yr1,Yr11,Yr2; pause Do iidx=1,Nidx Call rd_ncdf_2d (Fname,Year,Month,Nyr,Nmn,tmp, trim(Clim_nm(iidx))) j=0 Do iyr=Yr11,Yr2 Do i=1,Nyr if (Year(i)==iyr) then Do k=1,Nmon j=j+1; ClimIDX(j,iidx)=tmp(i,k) Enddo endif Enddo Enddo !pause Deallocate (Year,Month,tmp) !print "(12F7.3)",Y(iidx,:) !print*,Clim_nm(iidx);pause Enddo End Subroutine Read_ClimIDX !-------------------------------------------------------------------- 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)).OR.ABS(X(i))>1000) then Empty=.True.; Exit endif Enddo Empty_Ar1=Empty End Function Empty_Ar1 !--------------------------------------------------------------------------------- Subroutine Read_Dr_IDX (Fname) Character (Len=*) :: Fname Real, allocatable :: X3d(:,:,:),stnid(:),time(:),xlen(:) Integer :: Nt,Nle,it1 it1=Nmon+1 Call rd_ncdf_3d (Fname,stnid,time,xlen,Nst,Nt,Nle,X3d,"spi") !Allocate (SPI(Nst,Ntime,Nlen), SPEI(Nst,Ntime,Nlen)) Allocate (Dr_IDX(Nst,Ntime,Nlen,Ndr_idx)) ! (Nst,Ntime,Nlen,NDr_idx) Where (IsNaN(X3d)) X3d=Miss !SPI(1:Nst,1:Ntime,1:Nlen) = X3d(1:Nst,it1:Nt,1:Nlen) Dr_IDX(1:Nst,1:Ntime,1:Nlen,1) = X3d(1:Nst,it1:Nt,1:Nlen) Deallocate (stnid,time,xlen,X3d) Call rd_ncdf_3d (Fname,stnid,time,xlen,Nst,Nt,Nle,X3d,"spei") Where (IsNaN(X3d)) X3d=Miss !SPEI(1:Nst,1:Ntime,1:Nlen) = X3d(1:Nst,it1:Nt,1:Nlen) Dr_IDX(1:Nst,1:Ntime,1:Nlen,2) = X3d(1:Nst,it1:Nt,1:Nlen) Deallocate (stnid,time,xlen,X3d) Where (Dr_IDX(:,:,:,2)==Miss) Dr_IDX(:,:,:,2)=Dr_IDX(:,:,:,1) Call rd_ncdf_1d (Fname,stnid,Nst,xlon,"lon") Deallocate (stnid) Call rd_ncdf_1d (Fname,stnid,Nst,xlat,"lat") Deallocate (stnid) Call rd_ncdf_1d (Fname,stnid,Nst,Reg,"reg") Deallocate (stnid) Call rd_ncdf_1d_Char (Fname,stnid,Nst,stname,"stname") Deallocate (stnid) print*,"Calling Read_Dr_IDX Done" End Subroutine Read_Dr_IDX !----------------------------------------------------------------------------------! Subroutine Print_test Real, Allocatable, Dimension (:,:) :: A, B Integer :: i,j,k Allocate (A(Nmode,Nmode),B(Nmode,Nmode)) A = Matmul(Transpose(EOFs),EOFs); B = Matmul(Transpose(PCs),PCs) print*," --- Eof' x Eof --" Do i=1,Nmode print"(I3,15F8.3)",i,(A(i,j),j=1,Nmode) Enddo print*," --- PC' x PC --" Do i=1,Nmode print"(I3,15F8.3)",i,(B(i,j),j=1,Nmode) Enddo Deallocate (A,B) pause print*,"------" Do i=1,Nmode print*,Eigenvalues(i)/Sum(Eigenvalues)*100 Enddo print*," --- PCs ---" Do i=1,Nx print"(I3,15F8.3)",i,(PCs(i,j),j=1,Nmode) Enddo End Subroutine Print_test !------------------------------------------------------------------------------ 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" Include "Sub_GrADS.f90" END !----------------------------------------------------------------------------------! Include "/home/tanpv/Mylibs/Student_Prob.f90"