Use NetCDF Real, Parameter :: Miss=-99.0 Integer, Parameter :: Yr1=1960, Yr2=2015, Nmon=12 !--------------------------------------------------------------------------------- Real, Allocatable :: spei(:,:,:), prec(:,:,:), Mask(:,:) Real, Allocatable :: Lon(:), Lat(:), spei2d(:,:) Integer :: Nlon, Nlat, Ntime, LandPoints Integer :: ilon, ilat, it, ix Integer :: Nt,Nx,Ncount,Mcount,switch, Standard, PC_Stard,Nmode Integer :: i,j,K Real*8 :: Factor Real*8, Allocatable :: X(:,:) ! X(Nt,Nx) Real*8, Allocatable :: Y(:,:), Xbar(:), Std(:) Real*8, Allocatable :: EOFs(:,:), PCs(:,:) Real*8, Allocatable :: EValues(:) Character (Len=255) :: F_pca,F_eof !--------------------------------------------------------------------------------- ! For PCA: switch=1; Standard=1; PC_Stard=0; Nmode=15 !--------------------------------------------------------------------------------- Call Read_Data ! Return spei, prec, Mask, Lon, Lat, Nlon, Nlat, Ntime ! PCA with spei: Nt=Ntime; Nx=LandPoints print*,"Nt,Nx=",Nt,Nx; pause if (switch==1) then Mcount=Nt; Ncount=Nx else if (switch==0) then Mcount=Nx; Ncount=Nt endif Allocate (X(Nt,Nx),Y(Nt,Nx)) Allocate (Xbar(Nx),Std(Nx)) Allocate (EValues(Mcount),EOFs(Nt,Nmode),PCs(Nx,Nmode)) !ix=0 !Do ilat=1,Nlat ! Do ilon=1,Nlon ! if (Mask(ilon,ilat)/=Miss) then ! ix=ix+1; X(1:Nt,ix)=spei(ilon,ilat,1:Nt) ! endif ! Enddo !Enddo X(1:Nt,1:Nx) = spei2d(1:Nt,1:Nx) !print*,"X=",X(1,1:Nx); pause Call PCA_2019 (X,Y,Nt,Nx,switch,Standard,PC_Stard,Xbar,Std,EValues,EOFs,PCs,Nmode) print*,EValues(1:Nmode) !print*,PCs F_pca="../NC/pca_spei.nc" F_eof="../NC/eof_spei.nc" Call Write_EOF_PCA (F_pca,F_eof,"spei") CONTAINS !--------------------------------------------------------------------------------- Subroutine Write_EOF_PCA (F_pca,F_eof,varn) Character (Len=*) :: F_eof, F_pca, varn Character (Len=255) :: var_ln, t_unit Integer :: i,j,k,it Real :: tmp, total Real, Allocatable :: PCA(:,:,:), xmode(:), EOF1(:,:), X1(:), Xtmp(:) Real, Allocatable :: ECumulate(:), xtime(:) 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(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(varn) t_unit="months since 1976-01-01 " !if (iidx==1) then Call wr_ncdf_3d (F_pca,lon,lat,xmode,nlon,nlat,Nmode, PCA, & "pc", trim(t_unit),trim(var_ln),"-", "lon", "lat","mode", Miss) !else ! Call wr_ncdf_3d_apd (F_pca,PCA,trim(IDX_nm(iidx)),nlon,nlat, & ! Nmode, trim(var_ln), "-", "lon", "lat", "mode", Miss) !endif !------------------------------------------------------------------------- Xtmp(1:Nmode) = EValues(1:Nmode) Call wr_ncdf_1d_apd (F_pca, Xtmp, "eig", Nmode, & "Eigen values of "//trim(varn), "-", "mode", Miss) total=Sum(EValues(1:Mcount)) Allocate (ECumulate(Nmode)) !tmp=0.0 !Do i=1,Mcount ! tmp=tmp+EValues(i) ! if (I<=Nmode) then ! ECumulate(i)=0 ! Do j=1,I ! ECumulate(i)=ECumulate(i)+EValues(j) ! Enddo ! endif !Enddo Do i=1,Nmode ECumulate(i)=Sum(EValues(1:i))/total ! ECumulate(i)/tmp; Xtmp(i)=EValues(i)/total Enddo Call wr_ncdf_1d_apd (F_pca, Xtmp, "var",Nmode, & "Percent Contribution of "//trim(varn),"-", "mode", Miss ) Xtmp(1:Nmode) = ECumulate(1:Nmode) Call wr_ncdf_1d_apd (F_pca, Xtmp,"cum",Nmode, & "Cumuluted Eigenvalues of "//trim(varn), "-", "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 of "//trim(varn) !if (iidx==1) then !Call wr_ncdf_2d (F_eof,xmode,xtime, Nmode,Nt, EOF1, & ! trim(IDX_nm(iidx)), & ! trim(var_ln), "-", "mode", "mon", Miss) Call wr_ncdf_2d (F_eof,xtime,xmode, Nt,Nmode, EOF1, & "eof", trim(var_ln), "-", "mon","mode", Miss) !else !Call wr_ncdf_2d_apd (F_eof,EOF1, trim(IDX_nm(iidx)), Nmode,Nt, & ! trim(var_ln), "-", "mode", "mon", Miss) ! Call wr_ncdf_2d_apd (F_eof,EOF1, trim(IDX_nm(iidx)), Nt,Nmode, & ! trim(var_ln), "-", "mon", "mode", Miss) !endif print*,"Write_EOF_PCA: Done !!!" End Subroutine Write_EOF_PCA !--------------------------------------------------------------------------------- Subroutine Read_Data Character (Len=255) :: Fname Integer :: ilon,ilat,it,it1, Nt, ix Real, Allocatable :: X3d(:,:,:), Y3d(:,:,:), xtime(:), x1d(:) Fname="/work/projects/G09/SEACLID/GPCC_CRU_SEA/SPEI_SPI_GPCC_CRU.nc" Call rd_ncdf_3d (Fname,lon,lat,xtime,nlon,nlat,nt,X3d,"spei12") Deallocate (lon,lat,xtime) Call rd_ncdf_3d (Fname,lon,lat,xtime,nlon,nlat,nt,Y3d,"pr") Deallocate (lon,lat) Call rd_ncdf_2d (Fname,lon,lat,nlon,nlat,Mask,"land") Where (IsNaN(X3d)) X3d=Miss Where (IsNaN(Y3d)) Y3d=Miss !print*,Mask Ntime=Nt-Nmon; it1=Nmon+1 ! Bo nam dau tien !print*,"Nt,Ntime,it1: ",Nt,Ntime,it1; pause Allocate (spei(nlon,nlat,Ntime),prec(nlon,nlat,Ntime),x1d(Ntime)) spei=Miss; prec=miss; LandPoints=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)/=Miss) then x1d(1:Ntime)=X3d(ilon,ilat,it1:Nt) if (FullDara(X1d,Ntime,Miss)) then spei(ilon,ilat,1:Ntime)=X3d(ilon,ilat,it1:Nt) prec(ilon,ilat,1:Ntime)=Y3d(ilon,ilat,it1:Nt) else Mask(ilon,ilat)=Miss ! print*,ilon,ilat; print"(50F6.1)", x1d(Ntime-30:Ntime) endif endif if (Mask(ilon,ilat)>0) LandPoints=LandPoints+1 Enddo Enddo Allocate (spei2d(Ntime,LandPoints)) ix=0 Do ilat=1,Nlat Do ilon=1,Nlon if (Mask(ilon,ilat)>0) then ix=ix+1 spei2d(1:Ntime,ix)=spei(ilon,ilat,1:Ntime) endif Enddo Enddo End Subroutine Read_Data !------------------------------------------------------------------------------ Logical Function FullDara(X,N,Miss) Integer :: N, i, k Real :: X(N), Miss Logical :: Full Full=.True.; k=0 Do i=1,N !if (X(i)==Miss) then; Full=.False.; exit; endif if (X(i)/=Miss) k=k+1 Enddo !FullDara=Full if (k==N) then; FullDara=.True.; else; FullDara=.False.; endif End Function FullDara !--------------------------------------------------------------------------------- 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 "/home/tanpv/Progs/Drou_VN/EOF_New/Sub_PCA_2019.f90" !--------------------------------------------------------------------------------- END