Use netCDF Integer, Parameter :: Nmon=12, Yr1=1981, Yr2=2014, Nyrs=(Yr2-Yr1)+1 Real, Parameter :: Miss=-99.0, missing_value = 9.999e+20 Real, Parameter :: Dlon=0.5, Dlat=0.5, Lon1=60.0, Lon2=290.0, Lat1=-20.0, lat2=60.0 Real*8, Allocatable :: xdat(:,:), ydat(:,:), xtb(:), sx(:) Real*8, Allocatable :: EOFs(:,:), PCs(:,:) Real*8, Allocatable :: Eigenvalues(:), ECumulate(:) Integer :: Nt,Nx,Ncount,Mcount,switch, Standard, PC_Stard,Nmode Integer :: Ntime,Nst,Nlon,Nlat, ilon,ilat,ist Real, Allocatable :: TC(:,:), Corr(:,:), Prob(:,:) Real, Allocatable :: EOF1(:,:,:) Real, Allocatable :: sst(:,:,:,:) ! (Nlon,Nlat,Nyrs,Nmon) Real, Allocatable :: xlon(:),xlat(:),xyrs(:),xmon(:) ! Nlon,Nlat,Nyrs,Nmon Character (Len=250) :: Fname, pca_out, eof_out Character (Len=10) :: vname, tcvar, varf !------------------------------------------------------------------------------------------- Fname="/home/tanpv/Data/TCs/VES_Monthly_TC_NTC_Count_jp_new.nc" tcvar="tc6m" Call Read_TC (Fname,tcvar) !vname="vor" !"h850"! "vws" !"prmsl" ! "sst" vname="vor" !vname="vws" !vname="h850" !vname="prmsl" !vname="sst" Fname="/home/tanpv/Data/CFSR0.5/Vars/"//trim(vname)//".nc" pca_out="pca_"//trim(vname)//".nc" eof_out="eof_"//trim(vname)//".nc" varf=vname if (trim(varf)=="h850") then varf="hgt" Call Read_var_4d (Fname,varf) else Call Read_var (Fname,varf) endif Nx=Nlon*Nlat; Nt=Nyrs; switch=1; Standard=1; PC_Stard=0; Nmode=5 !--------------------------------------------------------------------------------- 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 (EOF1(Nt,Nmon,Nmode)) Do imon=1, Nmon ist=0 Do ilat=1,Nlat ! Firt Lon, last Lat Do ilon=1,Nlon ist=ist+1; xdat(1:nt,ist)=sst(ilon,ilat,1:nt,imon) Enddo Enddo Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, & Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode) !Call PCA_2019 (X,Y,Nt,Nx,switch, Standard, PC_Stard, & ! XBAR,STD,EValues,EOF,PCs,Nmode) Do imod=1,Nmode EOF1(1:Nt,imon,imod)=EOFs(1:Nt,imod) Enddo !print*,EOFs; pause Call Write_EOF_PCA (pca_out, eof_out, imon, vname, Yr1) Enddo Fname="Corr_"//trim(tcvar)//"_eof_"//trim(vname)//".nc" Call Calc_TQ (Fname, vname) !------------------------------------------------------------------------------------------- CONTAINS !------------------------------------------------------------------------------------------- Subroutine Calc_TQ (Fname, vname) Character (Len=*) :: Fname, vname Integer :: Lag1, lag,N,imtc,imvar,ytc1,ytc2,yva1,yva2 Integer :: iyr, imon,ilon,ilat Character (Len=50) :: var_nm, t_unit, var_ln, var_unit, & lon_nm, lat_nm, z_nm, t_nm Real :: Corre, Probe Real, Allocatable :: X(:), Y(:), xmon(:), xmode(:) Allocate (X(Nyrs),Y(Nyrs)) Allocate (Corr(1:Nmon,1:Nmode), Prob(1:Nmon,1:Nmode)) Do imon=1,Nmon ! TC mon X(:)=Miss; X(1:Nt)=TC(1:Nt,imon) !print*,X Do imod=1,Nmode Y(:)=Miss; Y(1:Nt)=EOF1(1:Nt,imon,imod) !print*,Y; pause Call Correlation (X, Y, Nyrs, Corre, Probe, Miss) !print*,Corr, Prob; pause Corr (imon,imod)=Corre Prob (imon,imod)=Probe Enddo Enddo Allocate (xmon(Nmon),xmode(Nmode)) Do imod=1,Nmode; xmode(imod)=imod; Enddo Do imon=1,Nmon; xmon(imon)=imon; Enddo lon_nm="mon"; lat_nm="mod"; var_unit="-"; var_nm="cor"; var_ln="Correlation Coefficients - TC-"//trim(vname) Call wr_ncdf_2d (Fname,xmon,xmode,Nmon,Nmode,Corr,var_nm, var_ln, var_unit, lon_nm, lat_nm, Miss) var_nm="prob"; var_ln="P-values of t-test" Call wr_ncdf_2d_apd (Fname,Prob,var_nm,Nmon,Nmode, var_ln, var_unit, lon_nm, lat_nm, Miss) End Subroutine Calc_TQ !------------------------------------------------------------------------------------------- Subroutine Read_TC (Fname,vname) Character (Len=*) :: Fname, vname Real, Allocatable :: xyear(:),xmon(:), xtc(:,:) Integer :: Ny, Nm, Nt, it, ix, iy, imon, iyr Call rd_ncdf_2d (Fname,xyear,xmon,Ny,Nm,xtc,trim(vname)) Allocate (TC(Nyrs,Nmon)) iy=0 Do iyr=Yr1,Yr2 Do it=1,Ny if (int(xyear(it))==iyr) then iy=iy+1; TC(iy,1:Nmon)=xtc(it,1:Nmon) endif Enddo !print"(I5,12F5.1)",iy,(TC(iy,1:Nmon)) Enddo End Subroutine Read_TC !------------------------------------------------------------------------------------------- Real Function Vorti (u1,v1,u2,v2,lonu1,latu1,lonu2,latu2,lonv1,latv1,lonv2,latv2) ! v2 ! !-----!-----! ! ! ! ! ! u1 !-----!-----! u2 ! ! ! ! ! !-----!-----! ! v1 Real :: u1,v1,u2,v2,lonu1,latu1,lonu2,latu2,lonv1,latv1,lonv2,latv2 Real :: du, dv, dx, dy du=u2-u1; dv=v2-v1 dx=Distance(lonu1,latu1,lonu2,latu2) dy=Distance(lonv1,latv1,lonv2,latv2) Vorti=du/dy-dv/dx End Function Vorti !------------------------------------------------------------------------------------------- Real Function Distance (lon1,lat1,lon2,lat2) Real :: lon1,lat1,lon2,lat2 Real, Parameter :: R=6371000, pi=3.1415926536, radi=pi/180.0 Real :: dlon, dlat, dlon2, dlat2, a, c, dis lon1=lon1*radi; lat1=lat1*radi; lon2=lon2*radi; lat2=lat2*radi dlon=lon2-lon1; dlat=lat2-lat1; dlon2=dlon/2; dlat2=dlat/2 a=sin(dlat2)*sin(dlat2) + cos(lat1)*cos(lat2)*sin(dlon2)*sin(dlon2) c=2*atan2(sqrt(a),sqrt(1-a)); dis=R*c Distance=dis End Function Distance !------------------------------------------------------------------------------------------- Subroutine Read_uv (Fname,vname) Character (Len=*) :: Fname,vname Integer :: Ntime, Nx, Ny, Nlev, it, ix, iy, imon, iyr Real, Allocatable :: xtmp(:,:,:), u(:,:,:,:), v(:,:,:,:) Real, Allocatable :: lon(:),lat(:), lev(:), xt(:) Call rd_ncdf_4d (Fname,lon,lat,lev,xt,nx,ny,nlev,ntime,u,"u") !Call rd_ncdf_3d (Fname,lon,lat,xt,nx,ny,ntime, X, trim(vname)) !"sst") Nlon=Int((Lon2-Lon1)/Dlon)+1; Nlat=Int((lat2-Lat1)/Dlat)+1 ix=Int(Lon1/Dlon)+1; iy=Int((Lat1+90)/Dlat)+1 Allocate (xtmp(Nlon,Nlat,Ntime)) Allocate (xlon(Nlon),xlat(Nlat)) xtmp(1:Nlon,1:Nlat,1:Ntime) = u(ix:ix+Nlon,iy:iy+Nlat,1,1:Ntime) & -u(ix:ix+Nlon,iy:iy+Nlat,2,1:Ntime) xlon(1:Nlon) = lon(ix:ix+Nlon); xlat(1:Nlat) = lat(iy:iy+Nlat) Allocate (sst(Nlon,Nlat,Nyrs,Nmon),xyrs(Nyrs),xmon(Nmon)) Do iyr=Yr1,Yr2; xyrs(iyr-Yr1+1)=iyr; Enddo Do imon=1,Nmon; xmon(imon)=imon; Enddo it=0 Do iyr=1,Nyrs Do imon=1,Nmon it=it+1 sst(1:Nlon,1:Nlat,iyr,imon)=xtmp(1:Nlon,1:Nlat,it) Enddo Enddo Where (sst==missing_value) sst=Miss print*,nx,ny,ntime,Nlon,Nlat,ix,iy print*,xlon(1),xlon(Nlon),xlat(1),xlat(Nlat) print*,ix,ix+Nlon,iy,iy+Nlat End Subroutine Read_uv !------------------------------------------------------------------------------------------- Subroutine Read_var (Fname,vname) Character (Len=*) :: Fname,vname Integer :: Ntime, Nx, Ny, it, ix, iy, imon, iyr Real, Allocatable :: xtmp(:,:,:), X(:,:,:),lon(:),lat(:),xt(:) Call rd_ncdf_3d (Fname,lon,lat,xt,nx,ny,ntime, X, trim(vname)) !"sst") Nlon=Int((Lon2-Lon1)/Dlon)+1; Nlat=Int((lat2-Lat1)/Dlat)+1 ix=Int(Lon1/Dlon)+1; iy=Int((Lat1+90)/Dlat)+1 Allocate (xtmp(Nlon,Nlat,Ntime)) Allocate (xlon(Nlon),xlat(Nlat)) xtmp(1:Nlon,1:Nlat,1:Ntime) = X(ix:ix+Nlon,iy:iy+Nlat,1:Ntime) xlon(1:Nlon) = lon(ix:ix+Nlon); xlat(1:Nlat) = lat(iy:iy+Nlat) Allocate (sst(Nlon,Nlat,Nyrs,Nmon),xyrs(Nyrs),xmon(Nmon)) Do iyr=Yr1,Yr2; xyrs(iyr-Yr1+1)=iyr; Enddo Do imon=1,Nmon; xmon(imon)=imon; Enddo it=0 Do iyr=1,Nyrs Do imon=1,Nmon it=it+1 sst(1:Nlon,1:Nlat,iyr,imon)=xtmp(1:Nlon,1:Nlat,it) Enddo Enddo Where (sst==missing_value) sst=Miss print*,nx,ny,ntime,Nlon,Nlat,ix,iy print*,xlon(1),xlon(Nlon),xlat(1),xlat(Nlat) print*,ix,ix+Nlon,iy,iy+Nlat End Subroutine Read_var !------------------------------------------------------------------------------------------- Subroutine Read_var_4d (Fname,vname) Character (Len=*) :: Fname,vname Integer :: Ntime, Nx, Ny, Nz, it, ix, iy, imon, iyr Real, Allocatable :: xtmp(:,:,:), X(:,:,:,:),lon(:),lat(:),lev(:),xt(:) !Call rd_ncdf_3d (Fname,lon,lat,xt,nx,ny,ntime, X, trim(vname)) !"sst") Call rd_ncdf_4d (Fname,lon,lat,lev,xt,nx,ny,nz,ntime,X,trim(vname)) Nlon=Int((Lon2-Lon1)/Dlon)+1; Nlat=Int((lat2-Lat1)/Dlat)+1 ix=Int(Lon1/Dlon)+1; iy=Int((Lat1+90)/Dlat)+1 Allocate (xtmp(Nlon,Nlat,Ntime)) Allocate (xlon(Nlon),xlat(Nlat)) xtmp(1:Nlon,1:Nlat,1:Ntime) = X(ix:ix+Nlon,iy:iy+Nlat,1,1:Ntime) xlon(1:Nlon) = lon(ix:ix+Nlon); xlat(1:Nlat) = lat(iy:iy+Nlat) Allocate (sst(Nlon,Nlat,Nyrs,Nmon),xyrs(Nyrs),xmon(Nmon)) Do iyr=Yr1,Yr2; xyrs(iyr-Yr1+1)=iyr; Enddo Do imon=1,Nmon; xmon(imon)=imon; Enddo it=0 Do iyr=1,Nyrs Do imon=1,Nmon it=it+1 sst(1:Nlon,1:Nlat,iyr,imon)=xtmp(1:Nlon,1:Nlat,it) Enddo Enddo if (trim(vname)=="sst") then Where (sst/=missing_value) sst=sst-273.15 endif Where (sst==missing_value) sst=Miss print*,nx,ny,ntime,Nlon,Nlat,ix,iy print*,xlon(1),xlon(Nlon),xlat(1),xlat(Nlat) print*,ix,ix+Nlon,iy,iy+Nlat End Subroutine Read_var_4d !------------------------------------------------------------------------- Subroutine Write_EOF_PCA (F_pca, F_eof, iidx, vname, init) Character (Len=*) :: F_eof, F_pca, vname Character (Len=255) :: var_ln, t_unit, varnm Integer :: i,j,k,iidx,it, init ! init=yr1, iidx=imon 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: write (varnm,"(A,I2.2)") trim(vname),iidx 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(varnm) write(t_unit,"(A,I4.4,A)") "years since ",Yr1,"-01-01 " if (iidx==1) then Call wr_ncdf_3d (F_pca,xlon,xlat,xmode,nlon,nlat,Nmode, PCA, & trim(varnm), trim(t_unit),trim(var_ln), & "-", "lon", "lat","mode", Miss) else Call wr_ncdf_3d_apd (F_pca,PCA,trim(varnm),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(varnm)//"_eig", Nmode, & "Eigen values of "//trim(varnm), "-", "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(varnm)//"_var",Nmode, & "Percent Contribution of "//trim(varnm), & "-", "mode", Miss ) Xtmp(1:Nmode) = ECumulate(1:Nmode) Call wr_ncdf_1d_apd (F_pca, Xtmp, & trim(varnm)//"_cum",Nmode, & "Cumuluted Eigenvalues of "//trim(varnm), & "-", "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(varnm) 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,xyrs,xmode, Nt,Nmode, EOF1, & trim(varnm), & trim(var_ln), "-", "year","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(varnm), Nt,Nmode, & trim(var_ln), "-", "year", "mode", Miss) endif print*,"Write_EOF_PCA: Done !!!" End Subroutine Write_EOF_PCA !------------------------------------------------------------------------------------------- Include "/home/tanpv/Mylibs/Sub_PCA_2019.f90" 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" END Include "/home/tanpv/Mylibs/Student_Prob.f90" !-------------------------------------------------------------------------------------------