Program PCA_T_R_ENSO ! Chuong trinh thuc hien PCA cho T, R, ENSO Use netcdf Integer, Parameter :: Yr1=1980, Yr2=2014, Nyear=Yr2-Yr1+1, Nmon=12, Ntime=Nyear*Nmon Integer, Parameter :: Nst=131, N_ClimIDX=16 Real, Parameter :: Miss=-99.0 !------------------------------------------ Include "FTNT.INC" ! PARAMETER (NMAX = Ntime, MMAX = Nst) !------------------------------------------ Character (Len = 11), Dimension(Nst) :: TenTram, VarName ! Ten bien Character (Len = 11), Dimension(N_ClimIDX) :: CL_Name Character (Len = 11), Dimension(Ntime) :: Sample ! Su dung nhu ten thu tu mau Integer, Dimension(Nst) :: YY ! So thu tu nam trong tap mau Real, Allocatable :: Pr(:,:), T2(:,:), Lon(:), Lat(:), H(:), Reg(:) Real, Allocatable :: PDSI(:,:), T2m(:,:), Prec(:,:) Real, Dimension (N_ClimIDX, Ntime) :: ClimIDX Character (Len=255) :: Res_Dir, Fname, Varn, Var_ln, Var_unit, lon_n, lat_n !------------------------------------------ ! Monthly T2m & Prec real :: Pre_mon(Nst, Nmon), T2m_mon(Nst, Nmon), T_P(Nst,Nmon*2) !------------------------------------------ ! X(Ntime,Nvar) = PC(Ntime,K) x EOF'(K,Nvar) ! TaiTrong = EOF(Nvar,K) ! NhanTo = PC (Ntime,K) ! Real, Dimension (Nmax,Mmax) :: XX ! Data to be used in Factor Analysis of Prec and Temp Real, Dimension (Nmax,N_ClimIDX) :: XX1 ! Data to be used in Factor Analysis of Clim_IDXs real :: TriRieng(MMax), VTORieng(MMax,MMax) ! MOI COT LA MOT VECTO real :: Eigens(Mmax) ! Luu TriRieng tich luy real :: TaiTrong(MMax,MMax) ! MOI COT LA TAI TRONG CUA MOT NHAN TO REAL :: NhanTo(NMax,MMax) ! MOI COT LA MOT NHAN TO REAL :: TT_DAQUAY(MMAX,MMAX), NT_DAQUAY(NMAX,MMAX) Real :: SumEig Common /Factors/ TaiTrong, NhanTo, TT_DAQUAY, NT_DAQUAY ! TaiTrong(M,K), NhanTo(N,K) Common /CEigen/ NEig, TriRieng, VTORieng, SumEig Integer :: NEig_Max ! So tri rieng duoc giu lai, =0 neu tong tich luy tri rieng>=95% Integer :: ist !------------------------------------------------------------------------------------------ Real :: Xtq(Ntime,12) ! Luu cac thanh phan de tinh tuong quan: 12=3(PDSI)+3(T2m)+3(Pre)+3(Cli) Real :: MTTQ(12,12), P_value(12,12), Corr, Prob ! Ma tran TQ Real :: Xtmp(Ntime), Ytmp(Ntime) Real, Allocatable :: Corr_PDSI_EOF_Cli(:,:), Prob_PDSI_EOF_Cli(:,:) Real, Allocatable :: Corr_PDSI_CliIDX(:,:), Prob_PDSI_CliIDX(:,:) Real, Allocatable :: Xlons(:), Xlats(:) ! Test: Real :: TQ_Clim_EOF(N_ClimIDX,N_ClimIDX) ! Kiem tra y nghia cua cac Tai trong Integer :: i, j Character (Len=6), Parameter :: Xtq_Name(12)=(/"PDSI1","PDSI2","PDSI3","T2m_1","T2m_2","T2m_3", & "Pre_1","Pre_2","Pre_3","Cli_1","Cli_2","Cli_3"/) !------------------------------------------------------------------------------------------ Res_Dir="/home/tanpv/SEACLID_New/RES_Stn/" ! Reading PDSI Fname="/home/tanpv/SEACLID_New/RES/PDSI_Res/PDSI_Stn_No_Rain_Center.nc" Call Read_PDSI_Stn (Fname) ! Read PDSI, T2m, Prec(Nst,Ntime), and Reg(Nst) Deallocate (Reg) Fname="/home/tanpv/SEACLID_New/Data/OBS_Stn/R_T2m_1980-2014_mon_Inland_No_Rain_Center.nc" Call Read_Prec_T2m_nc (FName) ! Return: Pr, T2 (Nst,Ntime); Lon, Lat, H, Reg (Nst) Call Read_Clim_IDX ! Return: ClimIDX, CL_Name (N_ClimIDX) Fname=trim(Res_Dir)//"ENSO_Phases.nc" ! Write ONI as ENSO definition Call Write_ENSO (Fname) Call Calc_Fill_Up_Miss (Pr, Nst, Ntime, Miss) ! Fill Up Missing values by Long term mean values Call Calc_Fill_Up_Miss (T2, Nst, Ntime, Miss) Call Calc_Fill_Up_Miss (ClimIDX, N_ClimIDX, Ntime, Miss) Call Calc_T2_Pre_Monthly ! Return Pre_mon, T2m_mon(Nst, Nmon), T_P(Nst,Nmon*2) !------------------------------------------ XX=transpose(Pr); NEig_Max=3 ! Precipitation Call PTNT(XX,Ntime,Nst,NEig_Max) Do ist=1, NEig; Xtq(1:Ntime,2*NEig+ist)=NhanTo(1:Ntime,ist); Enddo ! Xtq(Ntime,2*NEig:3*NEig) TriRieng = TriRieng/SumEig; Call Eigen_Accumulated print*,Eigens(1:NEig) Fname=trim(Res_Dir)//"PCA_Prec.nc"; Varn="PCA"; Var_ln="PCA of Precipitation"; Var_unit="-"; lon_n="Mon"; lat_n="Eigen" Call Write_PCA_EOF (Fname, Varn, Var_ln, Var_unit, lon_n, lat_n, Miss) !(Fname, NhanTo, TriRieng, Ntime, NEig, Varn, Var_ln, Var_unit, Lon, Lat, H, Reg, lon_n, lat_n, Miss) Fname=trim(Res_Dir)//"EOF_Prec.nc"; Varn="EOF"; Var_ln="EOF of Precipitation"; Var_unit="-"; lon_n="Stn"; lat_n="Eigen" Call Write_PCA_EOF (Fname, Varn, Var_ln, Var_unit, lon_n, lat_n, Miss) !(Fname, TaiTrong, TriRieng, Nst, NEig, Varn, Var_ln, Var_unit, Lon, Lat, H, Reg, lon_n, lat_n, Miss) !-------------------------------------- XX=transpose(T2); NEig_Max=2 ! 2m-Tempereture Call PTNT(XX,Ntime,Nst,NEig_Max) Do ist=1, NEig; Xtq(1:Ntime,NEig+ist)=NhanTo(1:Ntime,ist); Enddo ! Xtq(Ntime,NEig:2*NEig) TriRieng = TriRieng/SumEig; Call Eigen_Accumulated Fname=trim(Res_Dir)//"PCA_T2m.nc"; Varn="PCA"; Var_ln="PCA of 2m-Tempereture"; Var_unit="-"; lon_n="Mon"; lat_n="Eigen" Call Write_PCA_EOF (Fname, Varn, Var_ln, Var_unit, lon_n, lat_n, Miss) !(Fname, NhanTo, TriRieng, Ntime, NEig, Varn, Var_ln, Var_unit, Lon, Lat, H, Reg, lon_n, lat_n, Miss) Fname=trim(Res_Dir)//"EOF_T2m.nc"; Varn="EOF"; Var_ln="EOF of 2m-Tempereture"; Var_unit="-"; lon_n="Stn"; lat_n="Eigen" Call Write_PCA_EOF (Fname, Varn, Var_ln, Var_unit, lon_n, lat_n, Miss) !(Fname, TaiTrong, TriRieng, Nst, NEig, Varn, Var_ln, Var_unit, Lon, Lat, H, Reg, lon_n, lat_n, Miss) !-------------------------------------- Do ist=1, Nst; XX(ist,1:Nmon)=T2m_Mon(ist,1:Nmon); Enddo ! Long term mean of monthly 2m-Tempereture NEig_Max=2 Call PTNT(XX,Nst,Nmon,NEig_Max) TriRieng = TriRieng/SumEig; Call Eigen_Accumulated Fname=trim(Res_Dir)//"PCA_T2m_Mon.nc"; Varn="PCA_Mon"; Var_ln="PCA of 2m-Tempereture"; Var_unit="-"; lon_n="Stn"; lat_n="Eigen" Call Write_PCA_EOF (Fname, Varn, Var_ln, Var_unit, lon_n, lat_n, Miss) !(Fname, NhanTo, TriRieng, Ntime, NEig, Varn, Var_ln, Var_unit, Lon, Lat, H, Reg, lon_n, lat_n, Miss) Fname=trim(Res_Dir)//"EOF_T2m_Mon.nc"; Varn="EOF_Mon"; Var_ln="EOF of 2m-Tempereture"; Var_unit="-"; lon_n="Mon"; lat_n="Eigen" Call Write_PCA_EOF (Fname, Varn, Var_ln, Var_unit, lon_n, lat_n, Miss) !(Fname, TaiTrong, TriRieng, Nst, NEig, Varn, Var_ln, Var_unit, Lon, Lat, H, Reg, lon_n, lat_n, Miss) !-------------------------------------- Do ist=1, Nst; XX(1:Nmon,ist)=T2m_Mon(ist,1:Nmon); Enddo ! Long term mean of monthly 2m-Tempereture PA2 NEig_Max=2 Call PTNT(XX,Nmon,Nst,NEig_Max) TriRieng = TriRieng/SumEig; Call Eigen_Accumulated Fname=trim(Res_Dir)//"PCA_T2m_Mon_PA2.nc"; Varn="PCA_PA2"; Var_ln="PCA of 2m-Tempereture"; Var_unit="-"; lon_n="Mon"; lat_n="Eigen" Call Write_PCA_EOF (Fname, Varn, Var_ln, Var_unit, lon_n, lat_n, Miss) !(Fname, NhanTo, TriRieng, Ntime, NEig, Varn, Var_ln, Var_unit, Lon, Lat, H, Reg, lon_n, lat_n, Miss) Fname=trim(Res_Dir)//"EOF_T2m_Mon_PA2.nc"; Varn="EOF_PA2"; Var_ln="EOF of 2m-Tempereture"; Var_unit="-"; lon_n="Stn"; lat_n="Eigen" Call Write_PCA_EOF (Fname, Varn, Var_ln, Var_unit, lon_n, lat_n, Miss) !(Fname, TaiTrong, TriRieng, Nst, NEig, Varn, Var_ln, Var_unit, Lon, Lat, H, Reg, lon_n, lat_n, Miss) !-------------------------------------- Do ist=1, Nst; XX(ist,1:Nmon)=Pre_Mon(ist,1:Nmon); Enddo ! Long term mean of monthly Precipitation NEig_Max=3 Call PTNT(XX,Nst,Nmon,NEig_Max) TriRieng = TriRieng/SumEig; Call Eigen_Accumulated Fname=trim(Res_Dir)//"PCA_Pre_Mon.nc"; Varn="PCA_Mon"; Var_ln="PCA of Precipitation"; Var_unit="-"; lon_n="Stn"; lat_n="Eigen" Call Write_PCA_EOF (Fname, Varn, Var_ln, Var_unit, lon_n, lat_n, Miss) !(Fname, NhanTo, TriRieng, Ntime, NEig, Varn, Var_ln, Var_unit, Lon, Lat, H, Reg, lon_n, lat_n, Miss) Fname=trim(Res_Dir)//"EOF_Pre_Mon.nc"; Varn="EOF_Mon"; Var_ln="EOF of Precipitation"; Var_unit="-"; lon_n="Mon"; lat_n="Eigen" Call Write_PCA_EOF (Fname, Varn, Var_ln, Var_unit, lon_n, lat_n, Miss) !(Fname, TaiTrong, TriRieng, Nst, NEig, Varn, Var_ln, Var_unit, Lon, Lat, H, Reg, lon_n, lat_n, Miss) !-------------------------------------- Do ist=1, Nst; XX(1:Nmon,ist)=Pre_Mon(ist,1:Nmon); Enddo ! Long term mean of monthly Precipitation PA2 NEig_Max=3 Call PTNT(XX,Nmon,Nst,NEig_Max) TriRieng = TriRieng/SumEig; Call Eigen_Accumulated Fname=trim(Res_Dir)//"PCA_Pre_Mon_PA2.nc"; Varn="PCA_PA2"; Var_ln="PCA of Precipitation"; Var_unit="-"; lon_n="Mon"; lat_n="Eigen" Call Write_PCA_EOF (Fname, Varn, Var_ln, Var_unit, lon_n, lat_n, Miss) !(Fname, NhanTo, TriRieng, Ntime, NEig, Varn, Var_ln, Var_unit, Lon, Lat, H, Reg, lon_n, lat_n, Miss) Fname=trim(Res_Dir)//"EOF_Pre_Mon_PA2.nc"; Varn="EOF_PA2"; Var_ln="EOF of Precipitation"; Var_unit="-"; lon_n="Stn"; lat_n="Eigen" Call Write_PCA_EOF (Fname, Varn, Var_ln, Var_unit, lon_n, lat_n, Miss) !(Fname, TaiTrong, TriRieng, Nst, NEig, Varn, Var_ln, Var_unit, Lon, Lat, H, Reg, lon_n, lat_n, Miss) !-------------------------------------- Do ist=1, Nst; XX(ist,1:Nmon*2)=T_P(ist,1:Nmon*2); Enddo ! 2m-Tempereture + Precipitation: T_P(Nst,Nmon*2) NEig_Max=2 Call PTNT(XX,Nst,Nmon*2,NEig_Max) TriRieng = TriRieng/SumEig; Call Eigen_Accumulated Fname=trim(Res_Dir)//"PCA_T_P.nc"; Varn="PCA_Mon"; Var_ln="PCA of 2m-Tempereture+Precipitation"; Var_unit="-"; lon_n="Stn"; lat_n="Eigen" Call Write_PCA_EOF (Fname, Varn, Var_ln, Var_unit, lon_n, lat_n, Miss) !(Fname, NhanTo, TriRieng, Ntime, NEig, Varn, Var_ln, Var_unit, Lon, Lat, H, Reg, lon_n, lat_n, Miss) Fname=trim(Res_Dir)//"EOF_T_P.nc"; Varn="EOF_Mon"; Var_ln="EOF of 2m-Tempereture+Precipitation"; Var_unit="-"; lon_n="Mon"; lat_n="Eigen" Call Write_PCA_EOF (Fname, Varn, Var_ln, Var_unit, lon_n, lat_n, Miss) !(Fname, TaiTrong, TriRieng, Nst, NEig, Varn, Var_ln, Var_unit, Lon, Lat, H, Reg, lon_n, lat_n, Miss) !-------------------------------------- XX=transpose(PDSI); NEig_Max=3 ! PDSI Call PTNT(XX,Ntime,Nst,NEig_Max) Do ist=1, NEig; Xtq(1:Ntime,ist)=NhanTo(1:Ntime,ist); Enddo ! Xtq(Ntime,1:NEig) TriRieng = TriRieng/SumEig; Call Eigen_Accumulated Fname=trim(Res_Dir)//"PCA_PDSI.nc"; Varn="PCA"; Var_ln="PCA of PDSI"; Var_unit="-"; lon_n="Mon"; lat_n="Eigen" Call Write_PCA_EOF (Fname, Varn, Var_ln, Var_unit, lon_n, lat_n, Miss) !(Fname, NhanTo, TriRieng, Ntime, NEig, Varn, Var_ln, Var_unit, Lon, Lat, H, Reg, lon_n, lat_n, Miss) Fname=trim(Res_Dir)//"EOF_PDSI.nc"; Varn="EOF"; Var_ln="EOF of PDSI"; Var_unit="-"; lon_n="Stn"; lat_n="Eigen" Call Write_PCA_EOF (Fname, Varn, Var_ln, Var_unit, lon_n, lat_n, Miss) !(Fname, TaiTrong, TriRieng, Nst, NEig, Varn, Var_ln, Var_unit, Lon, Lat, H, Reg, lon_n, lat_n, Miss) !-------------------------------------- XX1=transpose(ClimIDX); NEig_Max=3 ! Clim_IDX Do ist=1, N_ClimIDX; XX(:,ist) = XX1(:,ist); Enddo ! Transfer XX1 to XX Call PTNT(XX,Ntime,N_ClimIDX,NEig_Max) Do ist=1, NEig; Xtq(1:Ntime,3*NEig+ist)=NhanTo(1:Ntime,ist); Enddo ! Xtq(Ntime,3*NEig:4*NEig) TriRieng = TriRieng/SumEig; Call Eigen_Accumulated Fname=trim(Res_Dir)//"PCA_ClimIDX.nc"; Varn="PCA"; Var_ln="PCA of Clim_IDXs"; Var_unit="-"; lon_n="Mon"; lat_n="Eigen" Call Write_PCA_EOF (Fname, Varn, Var_ln, Var_unit, lon_n, lat_n, Miss) !(Fname, NhanTo, TriRieng, Ntime, NEig, Varn, Var_ln, Var_unit, Lon, Lat, H, Reg, lon_n, lat_n, Miss) !Do ist=1, N_ClimIDX; print "(I4,5F10.4)", ist, TaiTrong(ist,1:NEig); Enddo Fname=trim(Res_Dir)//"EOF_ClimIDX.nc"; Varn="EOF"; Var_ln="EOF of Clim_IDXs"; Var_unit="-"; lon_n="IDX"; lat_n="Eigen" Call Write_PCA_EOF (Fname, Varn, Var_ln, Var_unit, lon_n, lat_n, Miss) !(Fname, TaiTrong, TriRieng, N_ClimIDX, NEig, Varn, Var_ln, Var_unit, Lon, Lat, H, Reg, lon_n, lat_n, Miss) !-------------------------------------- Fname=trim(Res_Dir)//"Corr_PDSI_ClimIDXs.nc" Call Calc_Corr_PDSI_ClimIDX (Fname) !-------------------------------------- ! Test: Tuong quan giua ClimIDXs va EOF cua ClimIDXs Do ist=1, N_ClimIDX Xtmp(1:Ntime)=ClimIDX(ist,1:Ntime) Do j=1, NEig Ytmp(1:Ntime)=NhanTo(1:Ntime,j) Call Correlation (Xtmp, Ytmp, Ntime, Corr, Prob, Miss) TQ_Clim_EOF(ist,j)=Corr Enddo print "(A15,10F10.3)", CL_Name(ist), (TQ_Clim_EOF(ist,j), j=1, NEig) Enddo !pause !-------------------------------------- ! Tuong quan giua EOF_Cli va PDSI tung tram: ! _PDSI_EOF_Cli(Nst,NEig) Allocate (Corr_PDSI_EOF_Cli(Nst, NEig), Prob_PDSI_EOF_Cli(Nst, NEig)) Do ist=1, Nst Xtmp(1:Ntime)=PDSI(ist,1:Ntime) Do j=1, NEig Ytmp(1:Ntime)=NhanTo(1:Ntime,j) Call Correlation (Xtmp, Ytmp, Ntime, Corr, Prob, Miss) Corr_PDSI_EOF_Cli(ist,j)=Corr; Prob_PDSI_EOF_Cli(ist,j)=Prob Enddo Enddo Allocate (Xlons(Nst), Xlats(NEig)) Do ist=1, Nst; Xlons(ist)=ist; Enddo; Do j=1, NEig; XLats(j)=Eigens(j); Enddo Fname=trim(Res_Dir)//"Corr_PDSI_EOF_Cli.nc" Varn="Corr"; Var_ln="Correlation between PDSI at each Station and EOF of ClimIDXs"; Var_unit="-"; lon_n="Stn"; lat_n="Eigen" Call wr_ncdf_2d (trim(Fname), Xlons, Xlats, Nst, NEig, Corr_PDSI_EOF_Cli, trim(Varn), trim(Var_ln), trim(Var_unit), trim(lon_n), trim(lat_n), Miss) Varn="Prob"; Var_ln="P-values of Correlation between PDSI at each Station and EOF of ClimIDXs"; Var_unit="-"; lon_n="Stn"; lat_n="Eigen" Call wr_ncdf_2d_apd (trim(Fname), Prob_PDSI_EOF_Cli, trim(Varn), Nst, NEig, trim(Var_ln), trim(Var_unit), trim(lon_n), trim(lat_n), Miss) Deallocate(Xlons,Xlats,Corr_PDSI_EOF_Cli,Prob_PDSI_EOF_Cli) !-------------------------------------- !-------------------------------------- Do i=1, 12 Do j=1, 12 Xtmp(1:Ntime)=Xtq(1:Ntime,i); Ytmp(1:Ntime)=Xtq(1:Ntime,j) Call Correlation (Xtmp, Ytmp, Ntime, Corr, Prob, Miss) MTTQ(i,j) = Corr; P_Value(i,j) = Prob Enddo Enddo Open (3,file="He_So_TQ.txt", Status="Unknown") Write (3,"(A)") " He so TQ giua cac EOF cua PDSI va T2m, Pre & ClimIDXs" Write (3,"(13A7)") " ",(Xtq_Name(i), i=1,12) Do i=1,12 Write (3,"(A6,1x,12F7.3)") Xtq_Name(i), (MTTQ(i,j),j=1,12) Enddo Write (3,"(A)") " P-value cua cac he so TQ giua cac EOF cua PDSI va T2m, Pre & ClimIDXs" Write (3,"(13A7)") " ",(Xtq_Name(i), i=1,12) Do i=1,12 Write (3,"(A6,1x,12F7.3)") Xtq_Name(i), (P_value(i,j),j=1,12) Enddo Close(3) !print*, NEig, (TriRieng(1:NEig)); pause !print*, VTORieng !------------------------------------------------------------------------------------------ CONTAINS Include "SUB/Calc_T2_Pre_Monthly.f90" Include "SUB/Calc_Corr_PDSI_ClimIDX.f90" Include "SUB/Write_ENSO.f90" Include "SUB/Write_PCA_EOF.f90" Include "SUB/Eigen_Accumulated.f90" Include "SUB/Calc_Fill_Up_Miss.f90" Include "SUB/Read_Prec_T2m_nc.f90" Include "SUB/Read_Clim_IDX.f90" Include "SUB/Read_PDSI_Stn.f90" Include "SUB/Time_Ave.f90" Include "/home/tanpv/SEACLID_New/Libs/Sub_Main_Factor.f90" Include "/home/tanpv/SEACLID_New/Libs/Sub_Factor.f90" Include "../Libs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/SEACLID_New/Libs/Sub_RW_Station_nc.f90" Include "../Libs/sub_rd_ncdf_Tan.F90" !Include "../Libs/Student_Test.f90" !Include "../Libs/Special_Function.f90" Include "SUB/Calc_Correlation_Coefficient.f90" !Include "../Libs/Student_Prob.f90" !------------------------------------------------------------------------------------------ END Include "../Libs/Student_Test.f90" Include "../Libs/Special_Function.f90" !Include "../Libs/Student_Prob.f90" !------------------------------------------------------------------------------------------