! Chuong trinh tinh PCs, EOFs theo cach moi (2019-10-10) Use NetCDF Implicit None Real, Parameter :: Miss=-99.0 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 Character (Len=250) :: Fname, Fname1, MaskFname, pca_out, eof_out Character (Len=10) :: vname Real*8 :: Factor Namelist /Model_Info/ Fname, vname, pca_out, eof_out, & Yr1, Yr2, Nmon, & switch, standard, PC_Stard, Nmode Integer :: i,j,L ! Mo ta thong tin SL dau vao: !--------------------------------------------------------------------------------- Real, allocatable :: xlon(:), xlat(:), xtime(:) Integer :: Nlon, Nlat, Ntime Logical :: Status Real, allocatable :: TC_tr(:,:,:), TC_Track(:,:) Integer :: Nst, iidx !--------------------------------------------------------------------------------- Open (1,file="Namelist_PCA_VES.txt",Status="old") Read (1, Model_Info) Close (1) !--------------------------------------------------------------------------------- ! Reading data: Call Read_TC_track (Fname,vname) Nx = Nst !--------------------------------------------------------------------------------- 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) = TC_Track(1:Nt,1:Nx) ! Tinh cac EOFs, PCs Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, & Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode) Call Write_EOF_PCA (pca_out,eof_out) print*, "===========> Complete " !Deallocate(Xbar,Std,Rx,EValues,V,Fx,Ax,D,EVtor,Lamda,EOF,PCs,ECumulate) print*,"Program is normal terminating ..." CONTAINS !================================================================================= Subroutine Write_EOF_PCA (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 Real, Allocatable :: PCA(:,:,:), xmode(:), EOF1(:,:), X1(:), Xtmp(:) 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 k=k+1; PCA(i,j,1:Nmode) = PCs(k,1:Nmode) Enddo Enddo Do i=1, Nmode; xmode(i)=i; Enddo ! PCA's modes var_ln="Spatial distribution of PCAs of TC tracks" t_unit="years since 1951-01-01 " Call wr_ncdf_3d (F_pca,xlon,xlat,xmode,nlon,nlat,Nmode, PCA, & "pc", trim(t_unit),trim(var_ln), & "-", "lon", "lat","mode", Miss) !------------------------------------------------------------------------- Xtmp(1:Nmode) = Eigenvalues(1:Nmode) Call wr_ncdf_1d_apd (F_pca, Xtmp, "eig", Nmode, & "Eigen values of ", "-", "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, "var",Nmode, & "Percent Contribution ", "-", "mode", Miss) Xtmp(1:Nmode) = ECumulate(1:Nmode) Call wr_ncdf_1d_apd (F_pca, Xtmp, "cum",Nmode, & "Cumuluted Eigenvalues", "-", "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" Call wr_ncdf_2d (F_eof,xtime,xmode, Nt,Nmode, EOF1, "eof", & trim(var_ln), "-", "year","mode", Miss) print*,"Write_EOF_PCA: Done !!!" End Subroutine Write_EOF_PCA !--------------------------------------------------------------------------------- Subroutine Read_TC_track (Fname,vname) Integer :: ist, ilon, ilat, it, it1, im, iyr, Ntime, Nyear Integer, Parameter :: Mtruot=5, Md2=Mtruot/2 Character (Len=*) :: Fname,vname Real, Allocatable :: X0(:), X1(:,:,:), X2(:,:,:), X3(:,:,:) Real, Allocatable :: Xmean(:), lon1(:),lat1(:) Integer :: nlon1, nlat1 Real :: tmp !Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, TC_tr, trim(vname)) ! "tc") Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, X1, trim(vname)) ! For VES: nlon1=nlon-8 ! 103E-123E Allocate (X2(nlon1,nlat,nt),lon1(nlon1)) lon1(1:Nlon1)=xlon(1:Nlon1) Do ilon=1,nlon1 Do ilat=1,nlat X2(ilon,ilat,1:nt)=X1(ilon,ilat,1:nt) Enddo Enddo Nlon=Nlon1 Deallocate (X1, xlon); Allocate (X1(Nlon,Nlat,Nt),xlon(Nlon)) X1(:,:,:)=X2(:,:,:); xlon(:)=lon1(:); Deallocate (X2,lon1) ! Done for VES Nyear=Yr2-Yr1+1; Allocate (TC_tr(nlon, nlat, nyear),X0(Nyear)) Do iyr=Yr1, Yr2 Do it=1, Nt if (Int(xtime(it))==iyr) then TC_tr(1:Nlon,1:Nlat,iyr-Yr1+1)=X1(1:Nlon,1:Nlat,it) X0(iyr-Yr1+1)=xtime(it) endif Enddo Enddo Deallocate (X1,xtime); Nt=Nyear Nst=Nlon*Nlat; Ntime=Nt-Mtruot+1 Allocate (xtime(Ntime)); xtime(1:Ntime)=X0(Md2+1:Nt-Md2) Deallocate (X0) Allocate (TC_Track(Ntime,Nst), Xmean(Nt)) Allocate (X1(Nlon,Nlat,Nt), X2(Nlon,Nlat,Ntime)) ! Areal Sum of TCs: Do it=1, Nt Xmean(it)=0 Do ilon=1, Nlon Do ilat=1, Nlat Xmean(it)=Xmean(it)+TC_tr(ilon,ilat,it) Enddo Enddo Enddo ! Convert to percentages: Do ilon=1, Nlon Do ilat=1, Nlat Do it=1, Nt X1(ilon,ilat,it)=TC_tr(ilon,ilat,it)/Xmean(it)*100 Enddo Enddo Enddo ! Smooth by running mean: Do ilon=1, Nlon Do ilat=1, Nlat it1=0 Do it=Md2+1, Nt-Md2 tmp=0; it1=it1+1 Do im=it-Md2,it+md2 tmp=tmp+X1(ilon,ilat,im) Enddo X2(ilon,ilat,it1)=tmp/Real(Mtruot) Enddo Enddo Enddo ! Convert to 2D-array: ist=0 Do ilat=1,Nlat Do ilon=1, Nlon ist=ist+1 TC_Track(1:Ntime,ist)=X2(ilon,ilat,1:Ntime) !TC_tr(ilon,ilat,1:Nt) Enddo Enddo Nt=Ntime Deallocate (Xmean,X1,X2) End Subroutine Read_TC_track !--------------------------------------------------------------------------------- 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))) 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