! 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 :: Nidx=24 Real*8, Allocatable :: xdat(:,:), ydat(:,:), xtb(:), sx(:) Real*8, Allocatable :: EOFs(:,:), PCs(:,:) Real*8, Allocatable :: Eigenvalues(:), ECumulate(:) !--------------------------------------------------------------------------------- Integer :: Yr1, Yr2, Nmon, 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) :: FnSST, pca_sst, eof_sst, pca_prec, eof_prec Real :: Slon1, Slon2, Slat1, Slat2, Miss_sst ! sst domain Integer :: Yr_sta Real*8 :: Factor Namelist /Model_Info/ Fname, MaskFname, pca_out, eof_out, & pca_prec, eof_prec, Yr1, Yr2, Nmon, & switch, standard, PC_Stard, Nmode Namelist /Clim_info/ Fnclim, pca_clim, eof_clim, & switch1, standard1, PC_Stard1, Nmode1 Namelist /SST_info/ FnSST, Yr_sta, Miss_sst, Slon1, Slon2, Slat1, & Slat2, pca_sst, eof_sst, & switch2, standard2, PC_Stard2, Nmode2 Integer :: i,j,L ! Mo ta thong tin SL dau vao: !--------------------------------------------------------------------------------- Real, allocatable :: xlon(:), xlat(:), xtime(:) Integer :: Nlon, Nlat, Ntime, ilon, ilat Logical :: Status Real, allocatable :: PDSI(:,:,:), SPI12(:,:,:), SPEI12(:,:,:) Real, allocatable :: ClimIDX(:,:) Real, allocatable :: Mask(:,:), Mask_all(:,:) Real, allocatable :: Dr_IDX(:,:,:) !--------------------------------------------------------------------------------- ! 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 Integer, Parameter :: NDr_idx=3 ! PDSI, SPEI12, SPI12 Character (Len=5), Parameter :: IDX_nm(NDr_idx)= & (/"pdsi","spi12","spe12"/) Character (Len=255) :: Filename !--------------------------------------------------------------------------------- !Character (Len=15), Parameter :: Clim_nm(Nidx)=(/ & ! "nino12", "nino3", "nino34", "nino4", "oni", "tni", & ! 6 !"whwp", "ammsst", "CAR_ersst", "NTA_ersst", "tsa", & ! 5 ! "whwp", "ammsst", "dmi", "NTA_ersst", "tsa", & ! 5 Bo CAR_ersst lay DMI ! "tna", "amonus", & ! "pna", "censo", "soi", "reqsoi", "np", "repac_slpa", & ! 6 ! "rindo_slpa", "qbo30", "wp", "ao", "nao"/) Character (Len=15), Parameter :: Clim_nm(Nidx)=(/ & "nino12", "nino3", "nino34", "nino4", "oni", "tni", & ! 6 "whwp", "ammsst", "pdo", "NTA_ersst", "dmi", & ! 5 Thay CAR_ersst bang pdo, thay tsa bang dmi "tna", "amonus", & "pna", "censo", "soi", "reqsoi", "np", "repac_slpa", & ! 6 "rindo_slpa", "qbo30", "wp", "ao", "nao"/) !--------------------------------------------------------------------------------- Open (1,file="Namelist_PCA_New.txt",Status="old") Read (1, Model_Info) Read (1, Clim_info) !Read (1, SST_info) Close (1) Mov_Ave_Step = 13 ! Moving average !--------------------------------------------------------------------------------- !Call Read_Drough_IDX (Fname) Nyear = Yr2-Yr1+1 ! So nam co SL !Nt = Ntime; Nx = Nst !--------------------------------------------------------------------------------- ! PCA for Clim: !--------------------------------------------------------------------------------- Call Read_ClimIDX (Fnclim) Nt = Ntime; Nx = Nidx switch=switch1; standard=standard1; PC_Stard=PC_Stard1; Nmode=Nmode1 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) 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" CONTAINS !================================================================================= 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 !------------------------------------------------------------------------- Subroutine Write_PCA_for_Clustering (Fname) Character (Len=*) :: Fname Real, Allocatable :: PCA(:,:),xmode(:),xpoint(:) ! (Nx,Nmode) Integer :: i,j,k,iidx,it 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_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 Write_EOF_PCA (F_pca, F_eof, iidx) Character (Len=*) :: F_eof, F_pca Character (Len=255) :: var_ln, t_unit, IDX_name, Vname Integer :: i,j,k,iidx,it Real :: tmp Real, Allocatable :: PCA(:,:,:), xmode(:), EOF1(:,:), X1(:), Xtmp(:) if (iidx>=1.AND.iidx<=3) then IDX_name=IDX_nm(iidx); Vname=IDX_nm(iidx) else if (iidx==0) then IDX_name="Precipitation"; Vname="prec" else IDX_name="SST or other variable"; Vname="sst" endif 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 if (Mask_all(i,j)>0) then k=k+1; PCA(i,j,1:Nmode) = PCs(k,1:Nmode) endif Enddo Enddo Do i=1, Nmode; xmode(i)=i; Enddo ! PCA's modes var_ln="Spatial distribution of PCAs of "//trim(IDX_name) write (t_unit,"(A,I4.4,A)") "months since ",Yr1,"-01-01 " if (iidx==1.OR.iidx<1.OR.iidx>3) then Call wr_ncdf_3d (F_pca,xlon,xlat,xmode,nlon,nlat,Nmode, PCA, & trim(Vname), trim(t_unit),trim(var_ln), & "-", "lon", "lat","mode", Miss) else Call wr_ncdf_3d_apd (F_pca,PCA,trim(Vname),nlon,nlat, & Nmode, trim(var_ln), "-", "lon", "lat", "mode", Miss) endif !------------------------------------------------------------------------- Xtmp(1:Nmode) = Eigenvalues(1:Nmode) Call wr_ncdf_1d_apd (F_pca, Xtmp, & trim(Vname)//"_eig", Nmode, & "Eigen values of "//trim(IDX_name), "-", "mode", Miss) tmp=0.0 Do i=1,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, & trim(Vname)//"_var",Nmode, & "Percent Contribution of "//trim(IDX_name), & "-", "mode", Miss ) Xtmp(1:Nmode) = ECumulate(1:Nmode) Call wr_ncdf_1d_apd (F_pca, Xtmp, & trim(Vname)//"_cum",Nmode, & "Cumuluted Eigenvalues of "//trim(IDX_name), & "-", "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(1:Nt,i)=EOFs(1:Nt,i); Enddo var_ln="Timeseries of eofs of "//trim(IDX_name) if (iidx==1.OR.iidx<1.OR.iidx>3) then Call wr_ncdf_2d (F_eof,xtime,xmode, Nt,Nmode, EOF1, & trim(Vname), & trim(var_ln), "-", "mon","mode", Miss) else Call wr_ncdf_2d_apd (F_eof,EOF1, trim(Vname), Nt,Nmode, & trim(var_ln), "-", "mon", "mode", Miss) endif print*,"Write_EOF_PCA: Done !!!" End Subroutine Write_EOF_PCA !--------------------------------------------------------------------------------- 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 !--------------------------------------------------------------------------------- !----------------------------------------------------------------------------------! 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/Student_Prob.f90"