! Chuong trinh tinh EOF cua 14 ClimIDX cho cac Lag tu 3-12 thang Use NetCDF ! Su dung SL cua cac thang tu 1979 den dau 2010 de Du bao cho cac thang tu 1/1980 den 12/2010 (Training) ! Real, Parameter :: Miss=-99.0 Integer, Parameter :: Yr1=1980, Yr2=2020, Nyear=Yr2-Yr1+1, Nmon=12, Ntime=Nyear*Nmon Integer, Parameter :: Nlen=2 Integer, Parameter :: Len(Nlen) = (/6,12/) Integer, Parameter :: Nidx=16 Character (Len=15), Parameter :: Clim_nm(Nidx)=(/"nino3", "nino34", "nino4", "pna", & !5 "censo", "oni", "soi", & "reqsoi", "repac_slpa", & !3 "rindo_slpa", "mei", "dmi", "pdo", & !4 "wpac85","cpac85", "epac85"/) !2 Real, Allocatable :: ClimIDX(:,:) ! (Nidx,Ntime Real, Allocatable :: SPI(:,:,:), SPEI(:,:,:) ! (Nst,Ntime,Nlen) Character(len=20), Allocatable :: stname(:) ! Nst Real, Allocatable :: Lon(:),Lat(:),Reg(:) !-------------------------------------------------------------------------------------- ! Tinh cac EOF cua cac ClimIDX: Real*8, Allocatable :: xdat(:,:), ydat(:,:), xtb(:), sx(:) Real*8, Allocatable :: EOFs(:,:), PCs(:,:) Real*8, Allocatable :: Eigenvalues(:), ECumulate(:) Real, Allocatable :: Xtb_Clim(:), Sx_Clim(:), PCs_Clim(:,:), Eig_Clim(:) Real, Allocatable :: lamda(:,:), PC1(:,:) Real, Allocatable :: EOF_Clim(:,:) Integer :: Nt, Nx ! Kich thuoc thuc cua mang X Integer :: Mcount, Ncount, switch, Standard, PC_Stard, Nmode, K Real*8 :: Factor Character (Len=200) :: FnameTC, FnameIDX, Fname, Out_Dir, Predictand !-------------------------------------------------------------------------------------- Integer :: imon, iyr, iidx, Nst Real, Allocatable :: Corr(:,:),Prob(:,:) ! (Nidx,Nidx) - Correlation between Clim_IDXa Real, Allocatable :: Cor_spi(:,:,:),Prb_spi(:,:,:) ! (Nst,Nlen,Nidx) - Correlation between Clim index and Dr index Real, Allocatable :: Cor_spei(:,:,:),Prb_spei(:,:,:) ! (Nst,Nlen,Nidx) Real, Allocatable :: Cor_spi_eof(:,:,:),Prb_spi_eof(:,:,:) ! (Nst,Nlen,Nmode) Real, Allocatable :: Cor_spei_eof(:,:,:),Prb_spei_eof(:,:,:) ! !-------------------------------------------------------------------------------------- FnameTC = "/home/tanpv/Progs/4BSc/Res/spi_spei_stn_1980-2020.nc" FnameIDX = "/work/projects/G09/SEACLID/ClimIDX_2022/Clim_IDX_2022.nc" Out_Dir = "/home/tanpv/Progs/4BSc/4Grads/" print*, "Calling Read_ClimIDX" Call Read_ClimIDX (FnameIDX) ! Do j=1,Ntime ! print "(I4,20F8.3)",j,(ClimIDX(iidx,j),iidx=1,Nidx) ! Enddo ! Correlations between climate indices: Call Calc_Correlation_Clim_Clim print "(A11,20A8)"," ",(Clim_nm(i),i=1,Nidx) Do i=1,Nidx print "(A8,20F8.3)",Clim_nm(i),(Corr(i,j),j=1,Nidx) Enddo ! Read drought indices: Call Read_Dr_IDX (FnameTC) ! Correlations between drought indices and climate indices: Call Calc_Correlation_Dr_Clim Fname=trim(Out_Dir)//"Cor_spi_ClimIDX" Call write_Correlation_Dr_Clim_to_GrADS (Cor_spi,Prb_spi,Fname) Fname=trim(Out_Dir)//"Cor_spei_ClimIDX" Call write_Correlation_Dr_Clim_to_GrADS (Cor_spei,Prb_spei,Fname) print "(A15,20A9)"," ",(trim(Clim_nm(i)),i=1,Nidx) Do i=1,Nlen Do j=1,Nst print "(A15,20F9.3)",trim(stname(j)),(Cor_spi(j,i,k),k=1,Nidx) Enddo pause Enddo !-------------------------------------------------------------------------------------- ! PCA of Climate indices: Call PCA_of_ClimIDX print*,"PC Clim:" print "(A15,20(A11,I1.1))"," ",("eof_",i,i=1,Nmode) Do j=1,Nidx print "(A15,20F12.3)",trim(Clim_nm(j)),(PCs_Clim(j,k),k=1,Nmode) Enddo pause ! Correlations between drought indices and eofs: Call Calc_Correlation_Dr_Eof_of_Clim Fname=trim(Out_Dir)//"Cor_spi_eof_of_ClimIDX" Call write_Correlation_Dr_eof_of_Clim_to_GrADS (Cor_spi_eof,Prb_spi_eof,Fname) Fname=trim(Out_Dir)//"Cor_spei_eof_of_ClimIDX" Call write_Correlation_Dr_eof_of_Clim_to_GrADS (Cor_spei_eof,Prb_spei_eof,Fname) print "(A15,20(A8,I1.1))"," ",("eof_",i,i=1,Nmode) Do i=1,Nlen Do j=1,Nst print "(A15,20F9.3)",trim(stname(j)),(Cor_spi_eof(j,i,k),k=1,Nmode) Enddo pause Enddo !-------------------------------------------------------------------------------------- CONTAINS !-------------------------------------------------------------------------------------- Subroutine PCA_of_ClimIDX Integer :: it, imode,iidx, i Real*8 :: Factor switch=0; Standard=1; PC_Stard=0; Nmode=6 Nt=Ntime; Nx=Nidx 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)) Allocate (EOF_Clim(Nt,Nmode)) Allocate (Xtb_Clim(Nx), Sx_Clim(Nx), PCs_Clim(Nx,Nmode), Eig_Clim(Nmode)) Do it=1,Nt !=Ntime Do i=1,Nx ! =Nidx*Nlag xdat(it,i)=ClimIDX(i,it) Enddo Enddo Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, & Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode) print"(A,6(I2,F6.2),F6.2)","Eigen - mode ",(imod, Eigenvalues(imod)/sum(Eigenvalues)*100,imod=1,Nmode),sum(Eigenvalues(1:Nmode))/sum(Eigenvalues)*100. Do imod=1,Nmode Do it=1,Nt EOF_Clim(it,imod)=EOFs(it,imod) Enddo Do iidx=1,Nx PCs_Clim(iidx,imod)=PCs(iidx,imod) Enddo Eig_Clim(imod)=Eigenvalues(imod) Enddo Do iidx=1,Nx Xtb_Clim(iidx)=Xtb(iidx); Sx_Clim(iidx)=Sx(iidx) Enddo End Subroutine PCA_of_ClimIDX !------------------------------------------------------------------------- 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, intent (in) :: 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 Calc_Correlation_Dr_Eof_of_Clim Integer :: i,j,k, Ntime1,lag Real, Allocatable :: X(:), Y(:) Real :: cor, prb Allocate (Cor_spi_eof(Nst,Nlen,Nmode),Cor_spei_eof(Nst,Nlen,Nmode)) Allocate (Prb_spi_eof(Nst,Nlen,Nmode),Prb_spei_eof(Nst,Nlen,Nmode)) Ntime1=Ntime-Nmon ! Bo nam 1980 do drought index khong du SL Allocate (X(Ntime1), Y(Ntime1)) lag = 6 ! so thang tre cua SL clim index Do i=1,Nmode X(1:Ntime1) = EOF_Clim(Nmon+1-lag:Ntime-lag,i) ! ClimIDX(i,Nmon+1-lag:Ntime-lag) Do k=1,Nlen Do j=1,Nst Y(1:Ntime1) = SPI(j,Nmon+1:Ntime,k) Call Correlation (X,Y,Ntime1,Cor,Prb,Miss) Cor_spi_eof(j,k,i)=Cor; Prb_spi_eof(j,k,i)=Prb Y(1:Ntime1) = SPEI(j,Nmon+1:Ntime,k) Call Correlation (X,Y,Ntime1,Cor,Prb,Miss) Cor_spei_eof(j,k,i)=Cor; Prb_spei_eof(j,k,i)=Prb Enddo Enddo Enddo End Subroutine Calc_Correlation_Dr_Eof_of_Clim !-------------------------------------------------------------------------------------- Subroutine Calc_Correlation_Dr_Clim Integer :: i,j,k, Ntime1,lag Real, Allocatable :: X(:), Y(:) Real :: cor, prb Allocate (Cor_spi(Nst,Nlen,Nidx),Cor_spei(Nst,Nlen,Nidx)) Allocate (Prb_spi(Nst,Nlen,Nidx),Prb_spei(Nst,Nlen,Nidx)) Ntime1=Ntime-Nmon ! Bo nam 1980 do drought index khong du SL Allocate (X(Ntime1), Y(Ntime1)) lag = 6 ! so thang tre cua SL clim index Do i=1,Nidx X(1:Ntime1) = ClimIDX(i,Nmon+1-lag:Ntime-lag) Do k=1,Nlen Do j=1,Nst Y(1:Ntime1) = SPI(j,Nmon+1:Ntime,k) Call Correlation (X,Y,Ntime1,Cor,Prb,Miss) Cor_spi(j,k,i)=Cor; Prb_spi(j,k,i)=Prb Y(1:Ntime1) = SPEI(j,Nmon+1:Ntime,k) !print*,Y(1:Ntime1); pause Call Correlation (X,Y,Ntime1,Cor,Prb,Miss) Cor_spei(j,k,i)=Cor; Prb_spei(j,k,i)=Prb Enddo Enddo Enddo End Subroutine Calc_Correlation_Dr_Clim !-------------------------------------------------------------------------------------- Subroutine write_Correlation_Dr_eof_of_Clim_to_GrADS (Cor,Prb,Fname) Character (Len=*) :: Fname Real, Dimension (Nst,Nlen,Nmode) :: Cor,Prb Real :: X(Nst,Nidx,Nlen*2) Integer :: ist,idx,Nlen2 Nlen2=Nlen*2 Do ist=1,Nst Do idx=1,Nmode X(ist,idx,1) = Cor(ist,1,idx) X(ist,idx,2) = Cor(ist,2,idx) X(ist,idx,3) = Prb(ist,1,idx) X(ist,idx,4) = Prb(ist,2,idx) Enddo Enddo Call Write_GrADS_Stn_Multivar (X,Nst,Nmode,Nlen2,Lon,Lat,stname,"spi",Fname) End Subroutine write_Correlation_Dr_eof_of_Clim_to_GrADS !-------------------------------------------------------------------------------------- Subroutine write_Correlation_Dr_Clim_to_GrADS (Cor,Prb,Fname) Character (Len=*) :: Fname Real, Dimension (Nst,Nlen,Nidx) :: Cor,Prb Real :: X(Nst,Nidx,Nlen*2) Integer :: ist,idx,Nlen2 Nlen2=Nlen*2 Do ist=1,Nst Do idx=1,Nidx X(ist,idx,1) = Cor(ist,1,idx) X(ist,idx,2) = Cor(ist,2,idx) X(ist,idx,3) = Prb(ist,1,idx) X(ist,idx,4) = Prb(ist,2,idx) Enddo Enddo Call Write_GrADS_Stn_Multivar (X,Nst,Nidx,Nlen2,Lon,Lat,stname,"spi",Fname) End Subroutine write_Correlation_Dr_Clim_to_GrADS !-------------------------------------------------------------------------------------- Subroutine Calc_Correlation_Clim_Clim Integer :: i,j Real, Allocatable :: X(:), Y(:) Real :: cor, prb Allocate (Corr(Nidx,Nidx),Prob(Nidx,Nidx)); Corr=Miss; Prob=Miss Allocate (X(Ntime), Y(Ntime)) Do i=1,Nidx X(1:Ntime) = ClimIDX(i,1:Ntime) Do j=1,Nidx Y(1:Ntime) = ClimIDX(j,1:Ntime) Call Correlation (X,Y,Ntime,Cor,Prb,Miss) Corr(i,j)=Cor; Prob(i,j)=Prb Enddo Enddo End Subroutine Calc_Correlation_Clim_Clim !-------------------------------------------------------------------------------------- Subroutine Read_ClimIDX (Fname) Character (Len=*) :: Fname Real, Allocatable :: Year(:),Month(:), tmp(:,:) Integer :: Nyr, Nmn, iidx, iyr, i,j, k, imon !Yr11=Yr1-1 ! Lui lai 1 nam Allocate (ClimIDX(Nidx,Ntime)) Do iidx=1,Nidx Call rd_ncdf_2d (Fname,Year,Month,Nyr,Nmn,tmp, trim(Clim_nm(iidx))) print*,"Reading ",iidx," ",trim(Clim_nm(iidx)) j=0 Do iyr=Yr1,Yr2 Do i=1,Nyr if (Year(i)==iyr) then Do imon=1,Nmon j=j+1; ClimIDX(iidx,j)=tmp(i,imon) !print "(F5.0,12F8.3)",Year(i) ,ClimIDX(iidx,j) !print*,j Enddo endif Enddo Enddo Deallocate (Year,Month,tmp) Enddo End Subroutine Read_ClimIDX !-------------------------------------------------------------------------------------- Subroutine Read_Dr_IDX (Fname) Character (Len=*) :: Fname Real, allocatable :: stnid(:),time(:),xlen(:) Integer :: Nt,Nle Call rd_ncdf_3d (Fname,stnid,time,xlen,Nst,Nt,Nle,SPI,"spi") Where (IsNaN(SPI)) SPI=Miss Deallocate (stnid,time,xlen) Call rd_ncdf_3d (Fname,stnid,time,xlen,Nst,Nt,Nle,SPEI,"spei") Where (IsNaN(SPEI)) SPEI=Miss Deallocate (stnid,time,xlen) Call rd_ncdf_1d (Fname,stnid,Nst,Lon,"lon") Deallocate (stnid) Call rd_ncdf_1d (Fname,stnid,Nst,Lat,"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) End Subroutine Read_Dr_IDX !-------------------------------------------------------------------------------------- Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Calc_Correlation_Coefficient.f90" Include "/home/tanpv/Mylibs/Sub_PCA_2019.f90" Include "/home/tanpv/Mylibs/Write_GrADS_Stn_Multivar.f90" !-------------------------------------------------------------------------------------- END !-------------------------------------------------------------------------------------- Include "/home/tanpv/Mylibs/Student_Prob.f90" !--------------------------------------------------------------------------------------