Use NetCDF Implicit None Real, Parameter :: Miss=-99.0 integer, parameter :: nmodes = 14 integer :: nt, nx, iminnxnt, imaxnxnt, nx_real, ny_real double precision, allocatable :: data(:,:), eigenvalues(:), & eigenvectors(:,:), & princomp(:,:), variance(:), & cumvariance(:) double precision :: value, time double precision, allocatable :: spacked(:), & evals(:), evecs(:,:), & rlawork(:), tentpcs(:,:), & sqrootweights(:), lon(:), & lat(:) integer, allocatable :: ilawork(:), & ifail(:) integer :: i, j, k, idx integer :: icovcor !------------------------------------------------------------------------------------ Real, allocatable :: GPI(:,:,:), xlon(:), xlat(:), xtime(:) Real, Allocatable :: Xmean(:,:), SS(:,:) Integer :: Nlon, Nlat, Ntime Character (Len=100) :: Fname, Fname1 Logical :: Status !------------------------------------------------------------------------------------ ! Test case: Should be commented after done nt=29; nx_real=180; ny_real=60 Ntime=nt; Nlon=nx_real; Nlat=ny_real !------------------------------------------------------------------------------ !Fname="../PCA/gpi_entire_domain.nc" !Call Read_GPI (Fname) !------------------------------------------------------------------------------ ! Identify the sizes of all arrays: Nt = Ntime; nx_real=Nlon; ny_real=Nlat Nx = nx_real*ny_real iminnxnt=min(nx,nt); imaxnxnt=max(nx,nt) ! Allocate: Allocate (data(nx,nt),eigenvalues(nmodes),eigenvectors(nx,nmodes), & princomp(nt,nmodes), variance(nmodes),cumvariance(nmodes)) Allocate (spacked(iminnxnt*(iminnxnt+1)/2),evals(iminnxnt), & evecs(iminnxnt,nmodes),rlawork(8*iminnxnt), & tentpcs(imaxnxnt,nmodes), & sqrootweights(nx), lon(nx_real), lat(ny_real), & ilawork(5*iminnxnt), ifail(iminnxnt)) Allocate (SS(nx,nx)) !------------------------------------------------------------------------------ ! Test case: Should be commented after done Call Read_test_data ! print *, 'read data done ......' !------------------------------------------------------------------------------ ! Convert 3D(lon,lat,t) to 2D-array (x,t): k=0 Do j=1, Nlat Do i=1, Nlon ! k=k+1; Data(k,1:Nt)=GPI(i,j,1:Nt) Enddo Enddo !------------------------------------------------------------------------------ do i=1, nx sqrootweights(i) = 1.0 ! Set area weights to 1.0 enddo icovcor = 0 ! use covariance matrix ! icovcor = 1 ! use correlation matrix Call deof( data, nx, nt, nmodes, icovcor,eigenvalues, & eigenvectors, & princomp, variance, cumvariance, & iminnxnt, imaxnxnt, spacked, evals, evecs, rlawork, & ilawork, ifail, tentpcs, sqrootweights ) !------------------------------------------------------------------------------ ! Testing: print*," FT x F =" print "(14F8.3)", matmul (transpose(eigenvectors),eigenvectors); pause print*," AT x A = L ??" print "(14F8.1)", transpose(matmul (transpose(princomp),princomp)); pause ! print "(14F8.3)", matmul (princomp,transpose(eigenvectors))-transpose(data); pause ! print "(14F8.1)", matmul (data,transpose(data)); pause ! print*," XT x F - A = 0" ! ===> OK ! print "(14F8.1)", matmul (transpose(data),eigenvectors)-princomp; pause ! print "(14F8.1)", princomp; pause print*," FT x XT x X x F = L - Diag - hope:" ! print "(14F8.1)", (matmul (matmul (transpose(eigenvectors),matmul (data,transpose(data))), eigenvectors)); pause print*, "XT x X = SS ??"; SS=matmul (data,transpose(data)) print "(14F8.1)", SS(1:14,1:14); pause !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ ! Write results to files: ! Fname="pca_eigen.nc"; Fname1="eof_timeseries.nc" Fname="pca_eigen_test.nc"; Fname1="eof_timeseries_test.nc" Call Write_EOF_PCA (Fname,Fname1) CONTAINS !------------------------------------------------------------------------------------ Subroutine Read_GPI (Fname) Character (Len=*) :: Fname Character (Len=255) :: var_ln Integer :: ilon, ilat, it Call rd_ncdf_3d (Fname, xlon, xlat, xtime, nlon, nlat, ntime, & GPI, "gpi") Allocate (Xmean(nlon, nlat)) Do ilon=1, Nlon Do ilat=1, Nlat Xmean(ilon,ilat)=0 Do it=1, Ntime Xmean(ilon,ilat) = Xmean(ilon,ilat) + GPI(ilon,ilat,it) Enddo Xmean(ilon,ilat) = Xmean(ilon,ilat)/Real(Ntime) Do it=1, Ntime GPI(ilon,ilat,it)=GPI(ilon,ilat,it)-Xmean(ilon,ilat) Enddo Enddo Enddo End Subroutine Read_GPI !----------------------------------------------------------------------------------! Subroutine Write_EOF_PCA (F_pca, F_eof) Character (Len=*) :: F_eof, F_pca Character (Len=255) :: var_ln Integer :: i,j,k Real, Allocatable :: PCA(:,:,:), xmode(:), EOF(:,:), X(:) !------------------------------------------------------------------------- Allocate (PCA(Nlon,Nlat,nmodes),xmode(nmodes)) k=0 Do j=1, Nlat Do i=1, Nlon k=k+1; PCA(i,j,1:nmodes) = eigenvectors(k,1:nmodes) Enddo Enddo Do i=1, nmodes; xmode(i)=i; Enddo ! PCA's modes var_ln="Spatial distribution of PCAs of GPI" Call wr_ncdf_3d (F_pca,xlon,xlat,xmode,nlon,nlat,nmodes, PCA, & "pca", "hours since 2010-08-01 ",trim(var_ln), "-", & "lon", "lat","mode", Miss) Call wr_ncdf_2d_apd (F_pca, Xmean, "mean",nlon,nlat, & "gpi mean", "-", "lon", "lat", Miss) Deallocate (PCA) !------------------------------------------------------------------------- ! Allocate (EOF(nmodes,Nt)) ! Do i=1, nmodes; EOF(i,1:Nt)=princomp(1:Nt,i); Enddo ! Call wr_ncdf_2d (F_eof,xmode,xtime, nmodes,Nt, EOF, "eof", & ! "Timeseries of eofs", "-", "mode", "time", Miss) ! Deallocate (EOF) Allocate (X(nmodes)) X(1:nmodes) = eigenvalues(1:nmodes) Call wr_ncdf_1d_apd (F_pca, X, "eigen", nmodes, & "Eigen values", "-", "mode", Miss) X(1:nmodes) = variance(1:nmodes) Call wr_ncdf_1d_apd (F_pca, X, "vari",nmodes, & "Variances", "-", "mode", Miss ) X(1:nmodes) = cumvariance(1:nmodes) Call wr_ncdf_1d_apd (F_pca, X, "cumvar",nmodes, & "Cumuluted variances", "-", "mode", Miss ) Deallocate (X) !------------------------------------------------------------------------- Allocate (EOF(nmodes,Nt)) Do i=1, nmodes; EOF(i,1:Nt)=princomp(1:Nt,i); Enddo Call wr_ncdf_2d (F_eof,xmode,xtime, nmodes,Nt, EOF, "eof", & "Timeseries of eofs", "-", "mode", "time", Miss) Deallocate (EOF) End Subroutine Write_EOF_PCA !----------------------------------------------------------------------------------! Subroutine Read_test_data ! nt=29; nx_real=180; ny_real=60; nmodes=13 print *, 'reading data' open(10,file='pacific_sst_test_data.txt',form='formatted',status='old') Allocate (xtime(nt)) do k=1, nt idx = 0; xtime(k)=k do j=1, ny_real do i=1, nx_real read(10,*) lon(i), lat(j), time, value idx = idx + 1 if( value .lt. -100.0 ) then data(idx,k) = 0.0 else data(idx,k) = value endif enddo enddo enddo Allocate (xlon(nx_real),xlat(ny_real),Xmean(nx_real,ny_real)) xlon(:)=lon(:); xlat(:)=lat(:); Xmean=1 End Subroutine Read_test_data !----------------------------------------------------------------------------------! Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" !------------------------------------------------------------------------------------ END