! Chuong trinh tinh EOF cua 14 ClimIDX cho cac Lag tu 3-12 thang Use NetCDF Real, Parameter :: Miss=-99.0 Integer, Parameter :: Yr1=1979, Yr2=2010, Nyear=Yr2-Yr1+1, Nmon=12, NTC_Cha=4 Integer, Parameter :: Yf1=2010, Yf2=2017, Nyfc= Yf2-Yf1+1 !Character (Len=8), Parameter :: Clim_nm(Nidx)=(/"ONI","MEI","REQSOI","SOI","RINDO","Nino34", & ! "BEST","PDO","PNA","DMI","WHWP","TNI","WP","QBO"/) !Character (Len=15), Parameter :: Clim_nm(Nidx)=(/"nino12", "nino3", "nino34", "nino4", "oni", "tni", & ! "mei", "pna", "censo", "ammsst", "CAR_ersst", "NTA_ersst", & ! "soi", "reqsoi", "whwp", "amonus", "np", "repac_slpa", & ! "rindo_slpa", "tsa", "tna","zwnd20", "wpac85", "epac85", & ! "cpac85", "qbo50", "qbo30", "wp", "ao", "nao"/) !Character (Len=5), Parameter :: TC_nm(NTC_Cha)=(/"tc","tc3m","tc6m","ntc","ntc3m","ntc6m"/) ! Nidx=24: !Integer, Parameter :: Nidx=24 !Character (Len=15), Parameter :: Clim_nm(Nidx)=(/"nino12", "nino3", "nino34", "nino4", "oni","tni",& ! "mei", "censo", "ammsst", "NTA_ersst", & ! "soi", "reqsoi", "whwp", "amonus", "repac_slpa", & ! "rindo_slpa", "tsa", "tna","zwnd20", "wpac85", "epac85", & ! "cpac85", "qbo30", "wp"/) ! Nidx=25: Integer, Parameter :: Nidx=25 Character (Len=15), Parameter :: Clim_nm(Nidx)=(/"nino12", "nino3", "nino34", "nino4", "oni", & !5 "pna", "censo", "ammsst", & !3 "soi", "whwp", "amonus", "tsa", "tna", "qbo30", "wp", "ao", & !8 "nao", "reqsoi", "repac_slpa", & !3 "rindo_slpa", "zwnd20", "wpac85", "epac85", & !4 "cpac85", "qbo50"/) !2 Character (Len=5), Parameter :: TC_nm(NTC_Cha)=(/"tc3m","tc6m","ntc3m","ntc6m"/) Real, Allocatable :: Clim_Dat(:,:,:) ! Nyear-1,Nlag*Nmon,Nmon (Nlag=10) Real, Allocatable :: Clim_Fcs(:,:,:), Y_Clim_Fcs(:,:) Real, Allocatable :: TC_Cha(:,:,:),TC_Fcs(:,:,:) Character (Len=200) :: FnameTC, FnameIDX, Fname, Out_Dir !-------------------------------------------------------------------------------------- ! Tinh cac EOF cua cac ClimIDX va tuong quan giua cac EOF voi TC_Char Real*8, Allocatable :: xdat(:,:), ydat(:,:), xtb(:), sx(:) Real*8, Allocatable :: EOFs(:,:), PCs(:,:) Real*8, Allocatable :: Eigenvalues(:), ECumulate(:) Integer :: Nt, Nx ! Kich thuoc thuc cua mang X Integer :: Mcount, Ncount, switch, Standard, PC_Stard, Nmode, K Real*8 :: Factor Real, Allocatable :: Xtb_Clim(:,:), Sx_Clim(:,:), PCs_Clim(:,:,:), Eig_Clim(:,:) Real, Allocatable :: lamda(:,:), PC1(:,:) Real, Allocatable :: EOF_Clim(:,:,:), ClimIDX(:,:,:) ! Loai bo CPAC850 va WPC850 vi SL ngan Real, Allocatable :: EOF_Fcs(:,:,:), EOF1(:,:), ClimFcs(:,:,:) Integer :: imon, iyr, iidx, imod, Nlag, Nt_Fcs Real, Allocatable :: Corr_all_lag(:,:,:),Prob_all_lag(:,:,:) ! Nmon,Nmode,NTC_Cha !-------------------------------------------------------------------------------------- ! So bien duoc chon boi Stepwise Regression Integer, Allocatable :: Mvar (:,:) ! (NTC_Cha,Nmon) Integer, Allocatable :: IDX_Var(:,:,:) ! (NTC_Cha,Nx,Nmon) Integer, Parameter :: Mvar_X=9 Real, Allocatable :: X_Clim(:,:,:,:), X_Fcs(:,:,:,:) Integer :: NMax, MMax, N, M Real :: PCT, ANS(11) Integer :: NSTEP(5), i,j, NN, KSTEP Integer, Allocatable :: IDX(:),L(:) Real, Allocatable :: X(:),XBAR(:),STD(:),D(:),B(:),RX(:),R(:),T(:) Integer, Allocatable :: NSTEPS(:,:),LS(:,:) Real, Allocatable :: BS(:,:), ANSS(:,:) Real, Allocatable :: XMEAN(:),DX(:), RR(:,:), YORG(:), YHQ(:) !-------------------------------------------------------------------------------------- FnameTC="/home/tanpv/Data/TCs/VES_Monthly_TC_NTC_Count_jp_new.nc" ! Contains: ntc, and tcday for 1,3,6 mon FnameIDX="/home/tanpv/Data/ClimIDX_2020/Clim_IDX_2020.nc" Call Read_TC (FnameTC) Call Read_ClimIDX (FnameIDX) Call Multi_lag ! Return: Clim_Dat switch=0; Standard=1; PC_Stard=0; Nmode=6 Nt=Nyear-1; Nlag=Nmon-2; Nx=Nidx*Nlag !-------------------------------------------------------------------------------------- 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 (EOF_Clim(Nt,Nmode,Nmon)) Allocate (Xtb_Clim(Nx,Nmon), Sx_Clim(Nx,Nmon), PCs_Clim(Nx,Nmode,Nmon), Eig_Clim(Nmode,Nmon)) !-------------------------------------------------------------------------------------- Do imon=1,Nmon print*,"=============imon==>",imon Do iyr=1,Nt !=Nyear Do iidx=1,Nx ! =Nidx xdat(iyr,iidx)=Clim_Dat(iyr,iidx,imon) if (xdat(iyr,iidx)==Miss) then; print*,"mon=",imon,xdat(iyr,iidx); endif Enddo !print"(14F10.3)",xdat(iyr,:);pause Enddo Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, & Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode) !print"(A,I2,8F8.2)","Eigen - Mon ",imon, Eigenvalues(1:Nmode), sum(Eigenvalues) !print"(A,I2,8F8.2)","Eigen - Mon ",imon, Eigenvalues(1:Nmode)/sum(Eigenvalues) print"(A,I2,8F8.2)","Eigen - Mon ",imon, sum(Eigenvalues(1:Nmode))/sum(Eigenvalues) !print"(14F8.3)",Ydat(:,:)/Factor; pause Do imod=1,Nmode Do iyr=1,Nt EOF_Clim(iyr,imod,imon)=EOFs(iyr,imod) Enddo Do iidx=1,Nx PCs_Clim(iidx,imod,imon)=PCs(iidx,imod) Enddo Eig_Clim(imod,imon)=Eigenvalues(imod) Enddo Do iidx=1,Nx Xtb_Clim(iidx,imon)=Xtb(iidx); Sx_Clim(iidx,imon)=Sx(iidx) Enddo Enddo ! Calculating EOF_Fcs: Nt_Fcs=Nyfc-1 ! Bo mot nam Allocate (Y_Clim_Fcs(Nt_Fcs,Nx),EOF_Fcs(Nt_Fcs,Nmode,Nmon),EOF1(Nt_Fcs,Nmode)) Allocate (lamda(Nmode,Nmode), PC1(Nx,Nmode)) lamda = 0.0 Do imon=1,Nmon Do iyr=1,Nt_Fcs Do iidx=1,Nx Y_Clim_Fcs(iyr,iidx)=(Clim_Fcs(iyr,iidx,imon)-Xtb_Clim(iidx,imon))/Sx_Clim(iidx,imon)/Factor Enddo !print"(14F8.3)",Y_Clim_Fcs(iyr,:) Enddo !pause PC1(1:Nx,1:Nmode)=PCs_Clim(1:Nx,1:Nmode,imon) !print"(6F8.3)",matmul(transpose(PC1),PC1); pause Do j=1,Nmode lamda(j,j) = 1./Eig_Clim(j,imon) ! 1/sqrt(eig) Enddo EOF1 = Matmul(Matmul(Y_Clim_Fcs,PC1),lamda) Do j=1,Nmode Do iyr=1,Nt_Fcs EOF_Fcs(iyr,j,imon) = EOF1(iyr,j) Enddo Enddo Enddo ! print test: Do imon=1,Nmon Do iyr=1,Nt print"(i4,10F8.3)",Yr1+iyr,EOF_Clim(iyr,:,imon) Enddo print* Do iyr=1,Nt_Fcs print"(i4,10F8.3)",Yf1+iyr,EOF_Fcs(iyr,:,imon) Enddo Enddo Call Create_File_for_Training ("Out_train_eof_cli/") ! Use EOFs Call Create_File_for_Forecast ("Out_fcs_eof_cli/") ! Use EOFs Call Var_Filter (Nt,Nx+1,Nt,Nx+1) Call Create_File_for_Training_cli ("Out_train_cli/") ! Use origin ClimIDX Call Create_File_for_forecast_cli ("Out_fcs_cli/") ! Use origin ClimIDX !-------------------------------------------------------------------------------------- CONTAINS !-------------------------------------------------------------------------------------- Subroutine Var_Filter (NMax, MMax, N, M) Integer :: NMax, MMax, N, M, KSTEP, i1, j1 Real :: PCT print*,"Calling Var_Filter ...", NMax, MMax, N, M Allocate (X(N*M),IDX(M),XMEAN(MMAX),DX(MMAX),RR(MMAX,MMAX),YORG(NMAX), YHQ(NMAX)) Allocate (NSTEPS(MMAX,5),LS(MMAX,MMAX),BS(MMAX,MMAX),ANSS(MMAX,11)) Allocate (Mvar(NTC_Cha,Nmon), IDX_Var(NTC_Cha,Mvar_X,Nmon)) Allocate (X_Clim(NTC_Cha,Nt,Mvar_X,Nmon), X_Fcs(NTC_Cha,Nt,Mvar_X,Nmon)) ! Mvar_X=9 Mvar=0; IDX_Var=Miss PCT=0.0000; IDX=0; IDX(1)=3; NS=1 Do itc=1,NTC_Cha Do imon=1,Nmon print*,"=============imon==>",imon X(1:Nt)=TC_Cha(itc,1:Nt,imon) !print*,"TC=",X(1:Nt); pause i=Nt Do iidx=1,Nx !iyr=1,Nt !=Nyear Do iyr=1,Nt !iidx=1,Nx ! =Nidx i=i+1; X(i)=Clim_Dat(iyr,iidx,imon) !print*,i,X(i) !xdat(iyr,iidx)=Clim_Dat(iyr,iidx,imon) Enddo Enddo !if (imon==3) then; print*,X; pause; endif Call HQ (NMax, MMax, N,M,NS,PCT,X,IDX, & XMEAN,DX, RR, YORG,YHQ, NSTEPS,LS,KSTEP,BS,ANSS) DO J=1,MMAX print*,j IF (NSTEPS(J,4).EQ.-999) THEN NN=J-1; EXIT ENDIF ENDDO ! ANSS(I,8),ANSS(I,6): Tuong ung la Error S, MulCoef R Mvar(itc,imon)=NN ! So bien duoc chon IDX_Var(itc,1:NN,imon)=LS(NN,1:NN) ! Chi so cac bien duoc chon print"(A,20I6)","1_imon,NN,ivar: ",imon,NN,IDX_Var(itc,1:NN,imon); print"(A,20I6)","LS:",imon,LS(NN,1:NN) print*,"Error, R:", ANSS(NN,8),ANSS(NN,6); !pause ! Sap xep lai thu tu chi so cac bien: Do i=1,NN-1 Do j=i+1,NN if (IDX_Var(itc,j,imon)Nidx) i1=Nidx endif !j1=mod(IDX_Var(itc,i,imon),10);if (j1==0) j1=10 lag=mod(IDX_Var(itc,i,imon),10); if (lag==0) lag=10 j1=imon-lag+10; if(j1>Nmon) j1=j1-Nmon !i1=IDX_Var(itc,i,imon)/10+1; if (i1>Nidx) i1=Nidx !j1=mod(IDX_Var(itc,i,imon),10);if (j1==0) j1=10 !j1=j1+2 write(EOFName(imon,i),"(A1,I2.2,A1,I2.2,A1,I2.2)") "X",i1,"_",j1,"_",imon Enddo Enddo Do imon=1, Nmon write (Fname,"(A,A,I2.2,A)") trim(Out_Dir)//trim(vnametc),"_",imon,".txt" write (Fname1,"(A,A,I2.2,A)") trim(Out_Dir)//"Ann_Out_"//trim(vnametc),"_",imon,".txt" write (Fname2,"(A,A,I2.2,A)") trim(Out_Dir)//"Ann_Save_"//trim(vnametc),"_",imon,".txt" !Write (33,"(A)") trim(Fname); Write (44,"(A)") trim(Fname1); Write (55,"(A)") trim(Fname2) Open (3,file=trim(Fname), Status="Unknown") write (3,"(A,I2.2,A)") "Total "//trim(vnametc)//" of month ", imon, " and Predictors" write (3,"(3I6)") Yr1, Nt, Mvar(itc,imon)+1 ! Mvar_X+1 write (3,"(20I3)") 3, (0,ivar=1,Mvar(itc,imon)) ! Mvar_X) !(IDX(i), i=1,Mvar) write (3,"(20A10)") trim(vnametc),(trim(EOFName(imon,i)),i=1,Mvar(itc,imon)) ! Mvar_X) Do iyr=1,Nt ! EOFall(LT,Nmon,Nmode,Nt) write (3,"(20F10.4)") TC(iyr,imon), EOF(iyr,1:Mvar(itc,imon),imon) ! EOF(iyr,1:Mvar_X,imon) Enddo Close (3) Enddo !Close (33); Close (44); Close (55) End Subroutine Create_file_for_SM_ClimIDX !-------------------------------------------------------------------------------------- Subroutine Create_File_for_Training (Out_Dir) Character (Len=*) :: Out_Dir Real :: TC(Nt,Nmon) Integer :: imon,iyr,imod, itc !Out_Dir="Out_train_eof_cli/" Do itc=1,NTC_Cha TC(1:Nt,1:Nmon)=TC_Cha(itc,1:Nt,1:Nmon) Call Create_file_for_Statistical_Models (Out_Dir,EOF_Clim,TC,Nt,trim(TC_nm(itc)),Yr1+1) Enddo End Subroutine Create_File_for_Training !-------------------------------------------------------------------------------------- Subroutine Create_File_for_Forecast (Out_Dir) Character (Len=*) :: Out_Dir Real :: TC(Nt_Fcs,Nmon) Integer :: imon,iyr,imod, itc !Out_Dir="Out_fcs_eof_cli/" Do itc=1,NTC_Cha TC(1:Nt_Fcs,1:Nmon)=TC_Fcs(itc,1:Nt_Fcs,1:Nmon) Call Create_file_for_Statistical_Models (Out_Dir,EOF_Fcs,TC,Nt_Fcs,trim(TC_nm(itc)),Yf1+1) Enddo End Subroutine Create_File_for_Forecast !-------------------------------------------------------------------------------------- Subroutine Create_file_for_Statistical_Models (Out_Dir,EOF,TC,Nt,vartc,Yr1) !(Lag, Kmon) !------------------------------------------------- !...Structure of the input data file as following (ASCII file): ! Row 1: Title of the file (CHARACTER) ! Row 2: Two positive integer number: N (Number of samples), Mbien (Number of variables) ! Row 3: Mbien values of IDX arrays (IDX(Mbien)) ! Row 4: Name of variables (CHARACTER*15) ! From row 5 until the end of file (N rows), each row contains Mbien values (REAL) Integer, Allocatable :: IDX(:) !=(/3,0,0,0,0,0,0,0,0,0/) ! 1 bien phu thuoc, 9 bien doc lap (3 months x 3 EOF) Integer :: ieigen, imon, iyr, isrc, i, Yr1, Nt Character (Len=7) :: EOFName(1:NMon,1:Nmode) Character (Len=255) :: Fname, Fname1, Fname2 Character (Len=10) :: vnametc Character (Len=*) :: Out_Dir, vartc Real :: EOF(Nt,Nmode,Nmon), TC(Nt,Nmon) vnametc=trim(adjustl(vartc)) Do ieigen=1, Nmode Do imon=1, NMon write(EOFName(imon,ieigen),"(A3,I1,A1,I2.2)") "EOF",ieigen,"_",imon Enddo Enddo Do imon=1, NMon write (Fname,"(A,A,I2.2,A)") trim(Out_Dir)//trim(vnametc),"_",imon,".txt" write (Fname1,"(A,A,I2.2,A)") trim(Out_Dir)//"Ann_Out_"//trim(vnametc),"_",imon,".txt" write (Fname2,"(A,A,I2.2,A)") trim(Out_Dir)//"Ann_Save_"//trim(vnametc),"_",imon,".txt" !Write (33,"(A)") trim(Fname); Write (44,"(A)") trim(Fname1); Write (55,"(A)") trim(Fname2) Open (3,file=trim(Fname), Status="Unknown") write (3,"(A,I2.2,A)") "Total "//trim(vnametc)//" of month ", imon, " and Predictors" write (3,"(3I6)") Yr1, Nt, Nmode+1 write (3,"(20I3)") 3, (0,ivar=1,Nmode) !(IDX(i), i=1,Mvar) write (3,"(20A10)") trim(vnametc),(trim(EOFName(imon,i)),i=1,Nmode) Do iyr=1,Nt ! EOFall(LT,Nmon,Nmode,Nt) write (3,"(20F10.4)") TC(iyr,imon), EOF(iyr,1:Nmode,imon) Enddo Close (3) Enddo !Close (33); Close (44); Close (55) End Subroutine Create_file_for_Statistical_Models !-------------------------------------------------------------------------------------- Subroutine Read_ClimIDX (Fname) Character (Len=*) :: Fname Real, Allocatable :: Year(:),Month(:), tmp(:,:) Integer :: Nyr, Nmn, iidx, iyr, i,j, k Allocate (ClimIDX(Nidx,Nyear,Nmon), ClimFcs(Nidx,Nyfc,Nmon)) Do iidx=1,Nidx Call rd_ncdf_2d (Fname,Year,Month,Nyr,Nmn,tmp, trim(Clim_nm(iidx))) j=0 Do iyr=Yr1,Yr2 Do i=1,Nyr if (Year(i)==iyr) then j=j+1; ClimIDX(iidx,j,:)=tmp(i,:) print "(F5.0,12F8.3)",Year(i),ClimIDX(iidx,j,:) endif Enddo Enddo ! For Independent Forecast: j=0 Do iyr=Yf1,Yf2 Do i=1,Nyr if (Year(i)==iyr) then j=j+1; ClimFcs(iidx,j,:)=tmp(i,:) print "(F5.0,12F8.3)",Year(i),ClimFcs(iidx,j,:) endif Enddo Enddo !print*,trim(Clim_nm(iidx));pause Deallocate (Year,Month,tmp) !print*,Clim_nm(iidx);pause Enddo End Subroutine Read_ClimIDX !-------------------------------------------------------------------------------------- Subroutine Multi_lag Integer :: imon, iidx, lag, imon1, imon11, Nlag, k Real, Allocatable :: X(:), Y(:) Nlag=Nmon-2 Allocate (Clim_Dat(Nyear-1,Nidx*Nlag,Nmon),Clim_Fcs(Nyfc-1,Nidx*Nlag,Nmon)) Allocate(X(Nyear-1), Y(Nyear-1)) Do imon=1, Nmon ! iidx=1, Nidx k=0 Do iidx=1, Nidx ! imon=1, Nmon Do lag=3,Nmon ! 0, Nmon -> Lay tu lag=3-lag=12 -> 10 lags imon1=imon-lag if (imon1>0) then X(1:Nyear-1)=ClimIDX(iidx,2:Nyear,imon1); !Y(1:Nyear-1)=TC_Cha(ist,2:Nyear,imon); Nxy=Nyear-1 else imon11 = 12 - ABS(imon1) X(1:Nyear-1)=ClimIDX(iidx,1:Nyear-1,imon11); !Y(1:Nyear-1)=TC_Cha(ist,2:Nyear,imon); Nxy=Nyear-1 endif k=k+1; Clim_Dat(1:Nyear-1,k,imon)=X(1:Nyear-1) !print*,"k=",Clim_Dat(1:Nyear-1,k,imon);pause Enddo Enddo Enddo Deallocate(X, Y) ! For Independent Forecast: Allocate(X(Nyfc-1), Y(Nyfc-1)) Do imon=1, Nmon ! iidx=1, Nidx k=0 Do iidx=1, Nidx ! imon=1, Nmon Do lag=3,Nmon ! 0, Nmon -> Lay tu lag=3-lag=12 -> 10 lags imon1=imon-lag if (imon1>0) then X(1:Nyfc-1)=ClimFcs(iidx,2:Nyfc,imon1); !Y(1:Nyear-1)=TC_Cha(ist,2:Nyear,imon); Nxy=Nyear-1 else imon11 = 12 - ABS(imon1) X(1:Nyfc-1)=ClimFcs(iidx,1:Nyfc-1,imon11); !Y(1:Nyear-1)=TC_Cha(ist,2:Nyear,imon); Nxy=Nyear-1 endif k=k+1; Clim_Fcs(1:Nyfc-1,k,imon)=X(1:Nyfc-1) !print*,"k=",Clim_Fcs(1:Nyfc-1,k,imon);pause Enddo Enddo Enddo Deallocate(X, Y) End Subroutine Multi_lag !-------------------------------------------------------------------------------------- Subroutine Read_TC (Fname) Character (Len=*) :: Fname Real, Allocatable :: Year(:),Month(:), tmp(:,:) Integer :: Nyr, Nmn, icha, iyr, i,j Allocate (TC_Cha(NTC_Cha,Nyear,Nmon),TC_Fcs(NTC_Cha,Nyfc,Nmon)) Do icha=1,NTC_Cha Call rd_ncdf_2d (Fname,Year,Month,Nyr,Nmn,tmp, trim(TC_nm(icha))) j=0 Do iyr=Yr1+1,Yr2 Do i=1,Nyr if (Year(i)==iyr) then j=j+1; TC_Cha(icha,j,:)=tmp(i,:) endif Enddo Enddo ! For Dependent Forecast: j=0 Do iyr=Yf1+1,Yf2 ! Bo nam dau tien vi trung voi Yr2 Do i=1,Nyr if (Year(i)==iyr) then j=j+1; TC_Fcs(icha,j,:)=tmp(i,:) endif Enddo Enddo Deallocate (Year,Month,tmp) Enddo End Subroutine Read_TC !-------------------------------------------------------------------------------------- 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/Mylibs/Calc_Correlation_Coefficient.f90" !Include "/home/tanpv/Mylibs/Sub_HQTB.f90" !-------------------------------------------------------------------------------------- END !-------------------------------------------------------------------------------------- Include "/home/tanpv/Mylibs/Sub_HQTB.f90" !--------------------------------------------------------------------------------------