Program FA_ENSO_Monthly ! 1) PCA cho cac ClimIDX dung lam nhan to du bao TC-count: ! 2) Create file for HQTB: ! + TC-count (imon) ! + EOFs (imon-3, imon-4, imon-5) of ClimIDX of Use NetCDF Integer, Parameter :: NMAX=1000, MMAX=20 Real, Parameter :: Miss=-99. ! Cac bien Output: Real :: TriRieng(MMax),VTORieng(MMax,MMax),TaiTrong(MMax,MMax),NhanTo(NMax,MMax),TT_DAQUAY(MMax,MMax), NT_DAQUAY(NMax,MMax) Real :: XMean(MMax), Stdev(MMax), Corr(MMax,MMax) Real :: XX(Nmax, Mmax) Integer :: YY(Nmax) !---------------------------------------------------------------------------- ! Origin Data: Real, Allocatable :: Xdat(:,:,:) !ClimIDX(:,:,:) Real, Allocatable :: TC(:,:),TC1(:,:) ! YTCs(:,:,:) !---------------------------------------------------------------------------- Real, Allocatable :: EOFs(:,:,:), PCs(:,:,:), Eig(:,:) Real, Allocatable :: Mean(:,:), Std(:,:) Character (Len = 15), Allocatable :: VarName(:), Sample(:), Yname(:) Integer, Allocatable :: mon_x(:), yrs_x(:) Character (Len=3), Allocatable :: MonName(:) Integer :: Nth0 ! So gia tri rieng duoc giu lai !---------------------------------------------------------------------------- Integer :: yr, iyr, ist, imon, K_Eigen Real :: Con Character (Len=4) :: Var Character (Len = 255) :: Fname, Out_Dir, Out_Dir_Fcs Integer :: NYr1, NYr2, Nmon, NtrMx, NEigenMax, NTC_src Integer :: Nyear, Ntime !---------------------------------------------------------------------------- Integer :: NLTime, Leadtime, NMonths, Nspl, Mvar ! So Leadtime, So thang du bao, So mode, Dung luong mau, So bien va so bien doc lap (Mvar1=Mvar-1) Integer, Allocatable :: Mon_fc(:) ! Thang duoc DB Character (Len=10) :: vname_tc ! Ten bien TC (phu thuoc) Character (Len=10), Allocatable :: vname(:,:) ! Ten cac bien truong Integer, Parameter :: Yr1=1982, Yr2=2009, Nyrs=(Yr2-Yr1+1), Mvarx=14, Nmode=5 Integer, Allocatable :: Nyears(:), Mvars(:), Mvar1(:), Yr_star(:) !Namelist /Lag_Kmon/ Lag, Kmon !Namelist /Data_info_Training/ NYr1, NYr2, Nmon, NtrMx, NEigenMax, NTC_src !Namelist /ClimIDX_Name/ VarName !Namelist /TC_Source/ Yname !Namelist /Mon_Name/ MonName !Namelist /Output_Dir/ Out_Dir, Out_Dir_Fcs !Namelist /PCA_inp_Out_Fname/ ClimFname, TCFname, PCA_out_Fname !---------------------------------------------------------------------------- Fname="/home/tanpv/Progs/TCs/FCS/Dat_out/TC_PCA_Leadtine_3.txt" Call Read_data (Fname) !Open (1, file="NameList_PCA.txt", Status="Old") !Read (1, Data_info_Training); print*, NYr1, NYr2, Nmon, NtrMx, NEigenMax,NTC_src ! Nyear = (Nyr2-Nyr1+1); Ntime = Nyear*Nmon; print*, Nyear, Ntime ! Allocate (VarName(NtrMx), Sample(Nmax), Yname(NTC_src)) ! Allocate (ClimIDX(1:Nyear,1:Nmon,1:NtrMx)) !Allocate (YTCs(1:Nyear,1:Nmon,1:NTC_src)) Allocate (EOFs(1:Nspl,1:Nmode,1:NMonths), PCs(1:Mvarx,1:Nmode,1:NMonths), Eig(1:Nmode,1:NMonths)) Allocate (Mean(1:Mvarx,1:NMonths), Std(1:Mvarx,1:NMonths)) Allocate (Yr_star(NMonths)) EOFs=Miss; PCs=Miss; Eig=Miss; Mean=Miss; Std=Miss Call Calc_PCA_EOF print*, "PCA Monthly... Done!" Out_Dir="/home/tanpv/Progs/TCs/FCS/Out/" Call Create_file_for_Statistical_Models Call Write_PCs_Eig_to_NC ("/home/tanpv/Progs/TCs/FCS/Dat_out/PCs_Fcs_traing.nc") !----------------------------------------------------------------------- CONTAINS !----------------------------------------------------------------------- Subroutine Write_PCs_Eig_to_NC (Fname) Character (Len=*) :: Fname Character (Len=4) :: PCname Real, Allocatable :: X(:,:,:), Y1(:,:), Y2(:,:), Nsample(:), Xlon(:), Xlat(:), Xtime(:) Real, Allocatable :: Z(:), xvar(:) Character (Len=255) :: var_ln, t_unit Integer :: Mvar Allocate (X(NMonths,Mvarx,Nmode), Y1(NMonths,Mvarx), Y2(NMonths,Mvarx),xvar(NMonths)) Allocate (Xlon(NMonths), Xtime(Nmode), Xlat(Mvarx), Z(NMonths),Nsample(NMonths)) t_unit="years since 1982-01-01" Do ilon=1,NMonths; Xlon(ilon)=Mon_fc(ilon); Enddo Do ilat=1,Mvarx; Xlat(ilat)=ilat; Enddo Do it =1,Nmode; Xtime(it)=it; Enddo var_ln="Cac thanh phan chinh" X=Miss; Y1=Miss; Y2=Miss Do imon=1,NMonths Mvar=Mvar1(imon) Do ivar=1,Mvar ! NEigenMax X (imon,ivar,1:Nmode) = PCs(ivar,1:Nmode,imon) Y1(imon,ivar) = Mean(ivar,imon) Y2(imon,ivar) = Std (ivar,imon) Enddo Nsample(imon) = Nyears(imon) Enddo !Call wr_ncdf_3d(trim(Fname),XLon,XLat,Xtime, Nmon,Ntr,NEigenMax, X, "PC", trim(t_unit), trim(var_ln),"-","Mon","ClimIDX","PCs",Miss) Call wr_ncdf_3d(trim(Fname),XLon,XLat,Xtime,NMonths,Mvarx,Nmode,X,"pc",trim(t_unit), trim(var_ln),"-","mon","var","PCs",Miss) var_ln="Trung binh thang cua cac bien" Call wr_ncdf_2d_apd (trim(Fname), Y1, "mean",NMonths,Mvarx, trim(var_ln), "-", "mon", "var", Miss) var_ln="Do lech chuan thang cua cac bien" Call wr_ncdf_2d_apd (trim(Fname), Y2, "std",NMonths,Mvarx, trim(var_ln), "-", "mon", "var", Miss) Do iei=1, Nmode write (var_ln,"(A,I1)") "Gia tri rieng pc",iei write (PCname,"(A,I1)") "eig",iei Z(1:NMonths) = Eig(iei,1:NMonths) Call wr_ncdf_1d_apd (trim(Fname), Z, trim(PCname), NMonths, trim(var_ln),"-", "mon", Miss) Enddo var_ln="Dung luong mau" Call wr_ncdf_1d_apd (trim(Fname), Nsample, "nspl", NMonths, trim(var_ln), "year", "mon", Miss) var_ln="So bien" xvar(:)=Mvar1(:) Call wr_ncdf_1d_apd (trim(Fname), xvar, "mvar", NMonths, trim(var_ln), "-", "mon", Miss) Deallocate (X,Y1,Y2,Z,Xlon,Xlat,Xtime,Nsample,xvar) End Subroutine Write_PCs_Eig_to_NC !----------------------------------------------------------------------------------! Subroutine Create_file_for_Statistical_Models !(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) !Real, Dimension (1:Nyear,1:Nmon,1:NTC_src) :: YTCs !Real :: EOFs(1:Nyear,1:NEigenMax,1:Nmon) 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, Mvar Character (Len=7) :: EOFName(1:NMonths,1:Nmode) Character (Len=255) :: Fname, Fname1, Fname2 Character (Len=10) :: vnametc vnametc=trim(adjustl(vname_tc)) ! Make Namelist file for TC_Traning_models: !Open (33,file="FnameInp.txt", Status="Unknown") !Open (44,file="FnameOut.txt", Status="Unknown") !Open (55,file="FnameSave.txt", Status="Unknown") !Do isrc=1, NTC_src Do ieigen=1, Nmode Do imon=1, NMonths write(EOFName(imon,ieigen),"(A3,I1,A1,I2.2)") "EOF",ieigen,"_",imon Enddo Enddo Do imon=1, NMonths !Mvar=Mvar1(imon) 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 TC count of month ", Mon_fc(imon), " and Predictors" write (3,"(3I6)") Yr_star(imon), Nyears(imon), Nmode+1 write (3,"(20I3)") 3, (0,ivar=1,Nmode) !(IDX(i), i=1,Mvar) !write (3,"(20A15)") "TC_"//MonName(imon),(EOFName(mon_x(k),:,isrc),k=1,Kmon) !write (3,"(20A10)") trim(vnametc),(trim(vname(imon,ivar)),ivar=1,Mvar) write (3,"(20A10)") trim(vnametc),(trim(EOFName(imon,i)),i=1,Nmode) Do iyr=1,Nyears(imon) ! Nspl write (3,"(20F10.4)") TC1(iyr,imon), EOFs(iyr,1:Nmode,imon) !write (*,"(20F15.3)") TC1(iyr,imon), EOFs(iyr,1:Mvar,imon) Enddo Close (3) Enddo !Enddo Close (33); Close (44); Close (55) End Subroutine Create_file_for_Statistical_Models !------------------------------------------------------------------------------------------- Subroutine Read_data (Fname) ! Read data of Climate Indices from NetCDF and store in 3D array of ClimIDX(1:Nyear,1:Nmon,1:Ntr) ! where Ntr=14 is number of Climate indices Integer :: iidx, iyr, imon, year, ivar, it, it0 Character (Len=*) :: Fname Real, Allocatable :: X(:) Real :: tmp Integer :: yyy, Mvar Allocate (X(Mvarx)) Open (1, file=Fname,status="Old") !read (1,*) NLTime Read (1,*) Leadtime ! Lead time cua DB hien tai read (1,*) NMonths ! So thang dc DB !print*,Leadtime,NMonths,Nyrs; pause Allocate (Mon_fc(NMonths)) ! Cac thang duoc du bao Allocate (Nyears(NMonths)) ! Dung luong mau thuc te (se duoc dung khi tinh PCs) Allocate (Mvars(NMonths),Mvar1(NMonths)) Allocate (Xdat(Nyrs,Mvarx,NMonths)) Allocate (TC(Nyrs,NMonths)) Allocate (vname(NMonths,Mvarx)) Do imon=1,NMonths read (1,*) Mon_fc(imon) read (1,*) Nspl, Mvars(imon) ! Dung luong mau, so bien Mvar1(imon)=Mvars(imon)-1 ! So bien doc lap Mvar=Mvar1(imon) !print*,Nspl,Mvar,Mvarx; pause !if (.not.Allocated(YY)) Allocate (YY(Nspl)) !if (.not.Allocated(TC)) Allocate (TC(Nspl,NMonths)) !if (.not.Allocated(Xdat)) Allocate (Xdat(Nspl,Mvar1,NMonths)) ! SL bien DL !if (.not.Allocated(vname)) Allocate (vname(Mvar1)) ! Ten bien doc lap read (1,"(20A10)") vname_tc, (vname(imon,ivar),ivar=1,Mvar) Do iyr=1,Nspl read (1,*) tmp, (X(it),it=1,Mvar) !print*,tmp,(X(it),it=1,Mvar) TC(iyr,imon)=tmp; Do ivar=1,Mvar; Xdat(iyr,ivar,imon)=X(ivar); Enddo !print*,(Xdat(iyr,ivar,imon),ivar=1,Mvar) YY(iyr)=Yr1+iyr-1 Enddo Enddo Allocate (Sample(Nyrs)) Allocate (TC1(Nyrs,NMonths)) ! (se duoc dung khi tinh PCs) End Subroutine Read_data !----------------------------------------------------------------------- Subroutine Calc_PCA_EOF !(Fname) ! Calculate the PCA and EOF and store in the 2D arrays of TaiTrong, NhanTo, respectively ! for each month of the year ! EOFs(1:Nyear,1:NEigenMax,1:Nmon) Integer :: iyr, i, j, Nt, Mvar Real :: xtmp(Mvarx) Do imon=1,NMonths Print*,"Tinh cho thang: ", Mon_fc(imon) !write (Fmt,'(I2.2)') Mon_fc(imon) Mvar=Mvar1(imon) CON = 1. !----------------- cac bien ------------ Nt=0 Loop_iyr: Do iyr=1,Nspl xtmp(1:Mvar)=Xdat(iyr,1:Mvar,imon) Do i=1,Mvar if (xtmp(i)==Miss) Cycle Loop_iyr Enddo Nt=Nt+1; if (Nt==1) Yr_star(imon)=Yr1+iyr-1 ! Nam bat dau co SL write(Sample(Nt),'(I4)') YY(iyr) !print*,imon,Nt,Yr_star XX(Nt,1:Mvar) = Xdat(iyr,1:Mvar,imon) !ClimIDX(iyr,imon,1:Ntr) TC1(Nt,imon)=TC(iyr,imon) ! De ghi ra cho DB Enddo Loop_iyr ! Nt=Nspl if no missing data Nyears(imon)=Nt ! De ghi vao file cho Du bao Nth0 = Nmode print*, Nspl, Nt, Mvar, Nmode, Nth0 Call PTNT(XX,NMAX,MMAX, Nt, Mvar, CON, Nmode, & XMean,Stdev,Corr,K_Eigen,TriRieng,VTORieng,TaiTrong,NhanTo,TT_DAQUAY,NT_DAQUAY) ! For other models use: !EOFs(1:,1:NEigenMax,imon) = NhanTo(1:Nyear,1:NEigenMax) EOFs(1:Nt,1:Nmode,imon) = NhanTo(1:Nt,1:Nmode) ! For Forecast Models: PCs (1:Mvar,1:Nmode,imon) = TaiTrong(1:Mvar,1:Nmode) Eig(1:Nmode,imon) = TriRieng(1:Nmode) Mean(1:Mvar,imon) = XMean(1:Mvar) Std (1:Mvar,imon) = Stdev(1:Mvar) !PCs (1:NtrMx,1:NEigenMax,imon) = TaiTrong(1:Ntr,1:NEigenMax) !Eig(1:NEigenMax,imon) = TriRieng(1:NEigenMax) !Mean(1:NtrMx,imon) = XMean(1:Ntr) !Std(1:NtrMx,imon) = Stdev(1:Ntr) Enddo ! k=1,Nmon End Subroutine Calc_PCA_EOF !----------------------------------------------------------------------- Subroutine Write_PCA_to_TXT (Fname,K_Eigen,N_Sample,M_Var) Integer :: K_Eigen,N_Sample,M_Var Character (Len=*) :: Fname Open(4,file=trim(Fname), status="Unknown") write(4,"(I3,A,16F8.3)") Nmode, " tri rieng: ",(TriRieng(i),i=1,Nmode) write(4,*) "Cac vector rieng:" Do j=1,Nmode write(4,"(I3,A,16F8.3)") j," : ", (VTORieng(i,j),i=1,Mvarx) Enddo write(4,*) "Ma tran Tai trong nhan to (PCA):" Do i=1, Mvarx write(4,"(A11,16F8.4)") VarName(i), (TaiTrong(i,j),j=1,Nmode) Enddo write(4,*) "Ma tran Tai trong nhan to da quay (PCA):" Do i=1, Mvarx write(4,"(A11,16F8.4)") VarName(i), (TT_DAQUAY(i,j),j=1,Nmode) Enddo write(4,*) "Ma tran nhan to (EOF):" Do i=1, N_Sample write(4,"(A11,16F8.4)") Sample(i), (NhanTo(i,j),j=1,Nmode) Enddo write(4,*) "Ma tran nhan to da quay (EOF):" Do i=1, N_Sample write(4,"(A11,16F8.4)") Sample(i), (NT_DAQUAY(i,j),j=1,Nmode) Enddo Close(4) End Subroutine Write_PCA_to_TXT !----------------------------------------------------------------------- Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" END ! Main Prog !----------------------------------------------------------------------- !Include "Old/Factor_PCA_EOF.f90" Include "/home/tanpv/Mylibs/Factor_PCA_EOF.f90"