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  :: ClimIDX(:,:,:)
Real, Allocatable :: 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 :: NTram, IO, yr, Ntr, 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 :: Kmon, Lag

Integer :: mon_y
Integer :: Nyrs, yrs_y
Character (Len=100) :: ClimFname, TCFname, PCA_out_Fname

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

!----------------------------------------------------------------------------

  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:Nyear,1:NEigenMax,1:Nmon), PCs(1:NtrMx,1:NEigenMax,1:Nmon), Eig(1:NEigenMax,1:Nmon))
    Allocate (Mean(1:NtrMx,1:Nmon), Std(1:NtrMx,1:Nmon))
    Allocate (MonName(Nmon))
  Read (1, Mon_Name)
  Read (1, ClimIDX_Name); Read (1, TC_Source)
  Read (1, Lag_Kmon)
  Read (1, Output_Dir)
  Read (1, PCA_inp_Out_Fname)
  Close (1)

  Ntr = NtrMx; !Out_Dir="Out/" ClimFname, TCFname, PCA_out_Fname

  Call Read_Clim_IDX(trim(ClimFname)) ; print*, "Call Read_Clim_IDX ... Done!"
  Call Read_TC_data (trim(TCFname)) ; print*, "Call Read_TC_data ... Done!"


  !Fname=trim(Out_Dir)//"PCA_EOF_Mon_"
  Call PCA_for_Monthly_Data !(Fname)

  print*, "PCA Monthly... Done!"

  ! Lag=3; Kmon=3
  !Do imon=1,Nmon
  !    mon_y=imon
  !    Call Predictor_Selection (mon_y, Lag, Kmon, Nyrs, yrs_y)
  
  !    print "(A)"," imon, mon_x1, mon_x2, mon_x3,  Nyrs, yrs_y, yrs_x1, yrs_x2, yrs_x3"
  !    print "(20I7)", imon,mon_x(1:Kmon), Nyrs, yrs_y, yrs_x(1:Kmon)
  !Enddo

  Call Create_file_for_Statistical_Models (Lag, Kmon)
  Call Make_Namelist_for_TC_Traing_models
  Call Write_PCs_Eig_to_NC (trim(PCA_out_Fname)) !"PCs_and_Eig_Values.nc")

  print*, "All_Timeseries... Done!"
!-----------------------------------------------------------------------

CONTAINS
!-----------------------------------------------------------------------
Subroutine Make_Namelist_for_TC_Traing_models
  Character (Len=100) :: Fname
  Integer :: ifile, Nfil
  
  Nfil = NTC_src*Nmon
  
  Open (3, file="NameList_TC_Training.txt", Status="Unknown")
! For ANN:
  Write (3,"(A)") "&INFO"
  Write (3,"(A)") "Nhidden       =     8        ! Number of hidden layers"
  Write (3,"(A)") "learning_rate =   0.5"
  Write (3,"(A)") "momentum      =     0"
  Write (3,"(A)") "training_mode =    .T.       ! For training"
  Write (3,"(A)") "Nepoc         = 15000        ! Iterative training times"
  Write (3,"(A)") "/"
  Write (3,"(A)") 
! For all:
  Write (3,"(A)") "&NumFile"
  Write (3,"(A,I4,A)") "Nfile = ", Nfil, ",                ! Number of data files to be calculated"   
  Write (3,"(A)") "/"
  Write (3,"(A)") 
  Write (3,"(A)") "&FileName"
  Write (3,"(A)") "  FnameInp ="
  Open (1, file="FnameInp.txt", Status="Old")
  Do ifile=1, Nfil
     Read (1,"(A)") Fname; Write (3,"(A)") "   '"//trim(Fname)//"'," 
  Enddo
  Close (1)
  Write (3,"(A)")
  Write (3,"(A)") "  FnameOut ="
  Open (1, file="FnameOut.txt", Status="Old")
  Do ifile=1, Nfil
     Read (1,"(A)") Fname; Write (3,"(A)") "   '"//trim(Fname)//"',"
  Enddo
  Close (1)
  Write (3,"(A)")
  Write (3,"(A)") "  FnameSave ="
  Open (1, file="FnameSave.txt", Status="Old")
    Do ifile=1, Nfil
     Read (1,"(A)") Fname; Write (3,"(A)") "   '"//trim(Fname)//"',"
  Enddo
  Close (1)
  Write (3,"(A)") "/"
  Close (3)
End Subroutine Make_Namelist_for_TC_Traing_models
!-----------------------------------------------------------------------
Subroutine Predictor_Selection (mon_y, Lag, Kmon, Nyrs, yrs_y) ! Chon 3 thang truoc thoi diem du bao (mon_y)
   Integer :: mon_y, yrs_y, Nyrs, Lag, Kmon, im
!   Integer :: mon_x(:), yrs_x(:)

   if (.not.Allocated(mon_x)) Allocate (mon_x(Kmon))
   if (.not.Allocated(yrs_x)) Allocate (yrs_x(Kmon))
   Lag_Km =  Lag + Kmon
   ista = mon_y - Lag_Km
   if (ista<=0) then
      m_sta = Nmon + ista
      Do im=1, Kmon 
         !mon_x(im)=m_sta+im-1
         if (m_sta+im-1 <= Nmon) then
            mon_x(im)=m_sta+im-1
            yrs_x(im)=Nyr1
         else
            mon_x(im)=mod(m_sta+im-1,Nmon) ! Nmon-(m_sta+im-1)
            yrs_x(im)=Nyr1+1
         endif  
      Enddo
      Nyrs = Nyear - 1
      yrs_y = Nyr1+1
   else
      m_sta = ista
      Do im=1, Kmon
         mon_x(im)=m_sta+im-1
         yrs_x(im)=Nyr1
         Nyrs = Nyear
         yrs_y = Nyr1
      Enddo
   endif
End Subroutine Predictor_Selection
!----------------------------------------------------------------------------------!
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(:)
  Character (Len=255) :: var_ln, t_unit

  Allocate (X(Nmon,Ntr,NEigenMax), Y1(Nmon,Ntr), Y2(Nmon,Ntr), Xlon(Nmon), Xtime(NEigenMax), Xlat(Ntr), Z(Nmon),Nsample(Nmon))
  t_unit="months since 1961-01-01"

  Do ilon=1,Nmon; Xlon(ilon)=ilon; Enddo
  Do ilat=1,Ntr; Xlat(ilat)=ilat; Enddo
  Do it=1, NEigenMax; Xtime(it)=it; Enddo
  var_ln="Cac thanh phan chinh"
  Do imon=1, Nmon
     Do itr=1, Ntr ! NEigenMax
        X(imon,itr,1:NEigenMax) = PCs(itr,1:NEigenMax,imon)
        Y1(imon,itr)            = Mean(itr,imon)
        Y2(imon,itr)            = Std(itr,imon)
     Enddo
     Nsample(imon) = Nyear
  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)

  var_ln="Trung binh thang cua cac ClimIDXs"
  Call wr_ncdf_2d_apd (trim(Fname), Y1, "Mean", Nmon,Ntr, trim(var_ln), "-", "Mon", "ClimIDX", Miss)
  var_ln="Do lech chuan thang cua cac ClimIDXs"
  Call wr_ncdf_2d_apd (trim(Fname), Y2, "Std", Nmon,Ntr, trim(var_ln), "-", "Mon", "ClimIDX", Miss)

  Do iei=1, NEigenMax
     write (var_ln,"(A,I1)") "Gia tri rieng PC",iei
     write (PCname,"(A,I1)") "Eig",iei
     Z(1:Nmon) = Eig(iei,1:Nmon) 
     Call wr_ncdf_1d_apd (trim(Fname), Z, trim(PCname), Nmon, trim(var_ln),"-", "Mon", Miss)
  Enddo

  var_ln="Dung luong mau"
  Call wr_ncdf_1d_apd (trim(Fname), Nsample, "Nspl", Nmon, trim(var_ln), "year", "Mon", Miss) 

  Deallocate (X,Y1,Y2,Z,Xlon,Xlat,Xtime,Nsample)

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, Parameter :: Mvar=10 
Integer, Parameter :: IDX(Mvar)=(/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
Character (Len=7) :: EOFName(1:Nmon,1:NEigenMax,1:NTC_src)
!Character (Len=3) :: MonName(Nmon) !=(/"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"/)
Character (Len=255) :: Fname, Fname1, Fname2
Integer :: mon_y, mon0
Integer :: yrs_y, Nyrs, Lag, Kmon

! 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, NEigenMax
      Do imon=1, Nmon
         write(EOFName(imon,ieigen,isrc),"(A3,I1,A1,I2.2)") "EOF",ieigen,"_",imon
      Enddo
   Enddo
   print "(12A8)", EOFName(:,:,isrc)

   Do imon=1, Nmon
      write (Fname,"(A,A,I2.2,A)")  trim(Out_Dir)//trim(Yname(isrc)),"_",imon,".txt"
      write (Fname1,"(A,A,I2.2,A)") trim(Out_Dir)//"Ann_Out_"//trim(Yname(isrc)),"_",imon,".txt"
      write (Fname2,"(A,A,I2.2,A)") trim(Out_Dir)//"Ann_Save_"//trim(Yname(isrc)),"_",imon,".txt"
      Write (33,"(A)") trim(Fname); Write (44,"(A)") trim(Fname1); Write (55,"(A)") trim(Fname2)  
      Open (3,file=trim(Fname), Status="Unknown")
      mon0=mod(imon+1,Nmon); if (mon0==0) mon0=Nmon
     ! write (3,"(A,A,I4,A)") "Total TC count of 6 months in advance since ", MonName(mon0)//"/",Nyr1, " and Predictors"
      mon_y = imon
      Call Predictor_Selection (mon_y, Lag, Kmon, Nyrs, yrs_y)
      write (3,"(A,A,I4,A)") "Total TC count of 6 months in advance since ", MonName(mon0)//"/",yrs_y, " and Predictors"
      write (3,"(3I6)") yrs_y, Nyrs, Mvar
      write (3,"(20I3)") (IDX(i), i=1,Mvar)
      write (3,"(20A15)") "TC_"//MonName(imon),(EOFName(mon_x(k),:,isrc),k=1,Kmon) 
      Do iyr=1, Nyrs
         write (3,"(20F15.3)") YTCs(yrs_y-NYr1+iyr,imon,isrc), (EOFs(yrs_x(k)-NYr1+iyr,:,mon_x(k)),k=1,Kmon) 
      Enddo
      Close (3)
   Enddo  
Enddo
  Close (33); Close (44); Close (55)
End Subroutine Create_file_for_Statistical_Models
!----------------------------------------------------------------------------------!
Subroutine Read_TC_data (Fname)
! Read TC data (TC count of 6 months as Predictands) from NetCDF file and
! store in the 3D array of YTCs(1:Nyear,1:Nmon,1:NTC_src)
! where: Nyear: 1961-2017, Nmon=12, NTC_src: JMA_N15,JMA_S15,JWTC_N15,JWTC_S15
Integer :: iidx, iyr, imon, year, it, it0
Character (Len=*) :: Fname
Real, Allocatable :: X(:,:),Xlon(:),Xlat(:)
Integer :: Nlon, Nlat

!Fname="/home/tanpv/OBS/KC09/TC_6mon.nc"
Do iidx=1, NTC_src
   if (Allocated(X)) Deallocate (X)
   if (Allocated(Xlon)) Deallocate (Xlon)
   if (Allocated(Xlat)) Deallocate (Xlat)

   Call rd_ncdf_2d(trim(Fname),Xlon,Xlat,Nlon,Nlat, X, trim(Yname(iidx)))
   Do iyr=NYr1, NYr2
      Do it=1, Nlon
      if (Int(Xlon(it))==iyr) then
         it0=iyr-NYr1+1
         YTCs(it0,:,iidx)=X(it,:)
      endif
      Enddo
   Enddo
Enddo
End Subroutine Read_TC_data

!----------------------------------------------------------------------------------!
Subroutine Read_Clim_IDX (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, IO, it, it0
Character (Len=*) :: Fname
Real, Allocatable :: X(:,:),Xlon(:),Xlat(:)
Integer :: Nlon, Nlat

!Fname="/home/tanpv/OBS/KC09/Clim_IDX.nc"

Do iidx=1, Ntr
   if (Allocated(X)) Deallocate (X)
   if (Allocated(Xlon)) Deallocate (Xlon)
   if (Allocated(Xlat)) Deallocate (Xlat)

   Call rd_ncdf_2d(trim(Fname),Xlon,Xlat,Nlon,Nlat, X, trim(VarName(iidx)))
!   it=0
   Do iyr=NYr1, NYr2
      Do it=1, Nlon
      if (Int(Xlon(it))==iyr) then
         it0=iyr-NYr1+1 
         ClimIDX(it0,:,iidx)=X(it,:)
         YY(it0) = iyr
      endif
      Enddo
   Enddo
Enddo
End Subroutine Read_Clim_IDX
!-----------------------------------------------------------------------
Subroutine PCA_for_Monthly_Data !(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)
!Character (Len=*) :: Fname
!Character (Len=255) :: Fname1, Fmt
Integer :: iyr, i, j
  Do imon=1,Nmon ! Mmon=1,2,...,12 correspond to Jan, Feb,..., Dec
     Print*,"Tinh cho thang: ",imon
     write (Fmt,'(I2.2)') imon   !
!     Fname1=trim(Fname)//trim(Fmt)//".txt"
     CON = 1.
   !----------------- cac bien ------------
     Do iyr=1,Nyear
        write(Sample(iyr),'(I4)') YY(iyr)
        XX(iyr,1:Ntr) = ClimIDX(iyr,imon,1:Ntr)
     Enddo
     Nth0 = NEigenMax
     print*, Nyear, Ntr, Nth0
     Call PTNT(XX,NMAX,MMAX,Nyear,Ntr,CON, Nth0, & 
               XMean,Stdev,Corr,K_Eigen,TriRieng,VTORieng,TaiTrong,NhanTo,TT_DAQUAY,NT_DAQUAY)
   ! For other models use:
     EOFs(1:Nyear,1:NEigenMax,imon) = NhanTo(1:Nyear,1:NEigenMax)
   ! For Forecast Models:
     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 PCA_for_Monthly_Data 
!-----------------------------------------------------------------------
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)") K_Eigen, " tri rieng: ",(TriRieng(i),i=1,K_Eigen)
     write(4,*) "Cac vector rieng:"
     Do j=1,K_Eigen
        write(4,"(I3,A,16F8.3)") j," : ", (VTORieng(i,j),i=1,M_Var)
     Enddo
     write(4,*) "Ma tran Tai trong nhan to (PCA):"
     Do i=1, M_Var
        write(4,"(A11,16F8.4)") VarName(i), (TaiTrong(i,j),j=1,K_Eigen)
     Enddo 
     write(4,*) "Ma tran Tai trong nhan to da quay (PCA):"
     Do i=1, M_Var
        write(4,"(A11,16F8.4)") VarName(i), (TT_DAQUAY(i,j),j=1,K_Eigen)
     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,K_Eigen)
     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,K_Eigen)
     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"