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 :: Yobs(:,:,:), Lon(:), Lat(:), Reg(:)
!----------------------------------------------------------------------------

  Real, Allocatable :: EOFs(:,:), PCs(:,:), Eig(:)  !EOFs(Ntime,NEigenMax)
  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, Nst
  Integer :: Kmon, Lag, Leadtime

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

  Namelist /Lag_Kmon/ Leadtime, 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, Predictand_Fname, PCA_out_Fname

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

  Open (1, file="NameList_PCA.txt", Status="Old")
  Read (1, PCA_inp_Out_Fname)       ! Provide filenames of ClimFname, Predictand_Fname, PCA_out_Fname
  Read (1, Data_info_Training)      ! Provide NYr1, NYr2, Nmon, NtrMx, NEigenMax, NTC_src
    Allocate (MonName(Nmon))
    Allocate (VarName(NtrMx), Yname(NTC_src))
  Read (1, Mon_Name)                ! Provide MonName
  Read (1, ClimIDX_Name)            ! Provide Climate index variables (VarName)
  Read (1, TC_Source)               ! Provide list of Drought IDXs (Yname)
  Read (1, Lag_Kmon)
  Read (1, Output_Dir)
  Close (1)

  Call Read_Predictand_data (trim(Predictand_Fname)) ! Return Yobs(Nst,Ntime,NTC_src), Lon(Nst), Lat(Nst), Reg(Nst), Nst, Ntime 
  ! Notes: Nst  is number of Stations
  !        Ntr  is number of Clim IDX

    Nyear = (Nyr2-Nyr1+1);  !   Ntime = Nyear*Nmon; print*, Nyear, Ntime
    Ntr   = NtrMx;
    if (Ntime /= Nyear*Nmon) then
       Print*, "Wrong in the input parameter ...."
       Stop
    endif
 
    Allocate (Sample(Ntime))
    Allocate (ClimIDX(1:Nyear,1:Nmon,1:NtrMx))
    Allocate (EOFs(1:Ntime,1:NEigenMax), PCs(1:Ntr,1:NEigenMax), Eig(1:NEigenMax))
    Allocate (Mean(1:Ntr), Std(1:Ntr))
  !  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)

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

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

  Call Create_file_for_Statistical_Models (Lag, Kmon)
  Call Make_Namelist_for_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_Traing_models
  Character (Len=100) :: Fname
  Integer :: ifile, Nfil
  
  Nfil = NTC_src*Nst
  
  Open (3, file="NameList_Drought_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_Traing_models
!-----------------------------------------------------------------------
Subroutine Write_PCs_Eig_to_NC (Fname)
  Character (Len=*) :: Fname
  Character (Len=4) :: PCname
  Real, Allocatable :: X(:,:), Y1(:), Y2(:), Nsample(:), Xlon(:), Xlat(:)
  Real, Allocatable :: Z(:)
  Character (Len=255) :: var_ln, t_unit

  Allocate (X(Ntr,NEigenMax), Y1(Ntr), Y2(Ntr), Xlon(Ntr), Xlat(NEigenMax), Z(NEigenMax),Nsample(NEigenMax))
  write(t_unit,"(A,I4.4,A)") "months since ",NYr1,"-01-01"


  Do ilon=1,Ntr; Xlon(ilon)=ilon; Enddo        ! ClimIDXs
  Do ilat=1,NEigenMax; Xlat(ilat)=ilat; Enddo  ! Eigens
  var_ln="Cac thanh phan chinh"
  Do itr=1, Ntr
     X(itr,1:NEigenMax) = PCs(itr,1:NEigenMax)
     Y1(itr)            = Mean(itr)
     Y2(itr)            = Std(itr)
  Enddo
  Call wr_ncdf_2d (trim(Fname),Xlon,Xlat,Ntr,NEigenMax, X, "PC", trim(var_ln),"-","ClimIDX", "PCs",Miss)
  var_ln="Trung binh cua cac ClimIDXs"
  Call wr_ncdf_1d_apd (trim(Fname), Y1, "Mean", Ntr, trim(var_ln),"-","ClimIDX",Miss)
  var_ln="Do lech chuan cua cac ClimIDXs"
  Call wr_ncdf_1d_apd (trim(Fname), Y2, "Std", Ntr, trim(var_ln),"-","ClimIDX",Miss)

  Z(:) = Miss; Nsample(:) = Miss
  Z(1:NEigenMax) = Eig(1:NEigenMax); Nsample(1:NEigenMax) = Ntime
  var_ln="Cac gia tri rieng "
  Call wr_ncdf_1d_apd (trim(Fname), Z, "Eig", NEigenMax, trim(var_ln),"-","PCs", Miss)

  Nsample(:) = Ntime
  var_ln="Dung luong mau"
  Call wr_ncdf_1d_apd (trim(Fname), Nsample, "Nspl", NEigenMax, trim(var_ln),"Mon", "PCs", Miss)

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

End Subroutine Write_PCs_Eig_to_NC
!----------------------------------------------------------------------------------!
Subroutine Create_file_for_Statistical_Models (Lag, Kmon)
  Integer, Parameter :: Mvar=10
  Integer, Parameter :: IDX(Mvar)=(/3,0,0,0,0,0,0,0,0,0/) 
  Integer :: ieigen, imon, iyr, isrc, i
  Character (Len=7) :: EOFName(1:Mvar-1)  !(1:Kmon,1:NEigenMax,1:NTC_src)
  Character (Len=255) :: Fname0, Fname, Fname1, Fname2
  Integer :: Lag, Kmon, Nspl, ista
  Real, Allocatable :: X(:), Y(:)
  Integer :: mon_y, mon_x1, mon_x2, mon_x3, mon0
  Integer :: yrs_y, Nyrs, yrs_x1, yrs_x2, yrs_x3

  ! Make Namelist file for Traning_models:
   Open (33,file="FnameInp.txt", Status="Unknown")
   Open (44,file="FnameOut.txt", Status="Unknown")
   Open (55,file="FnameSave.txt", Status="Unknown")

   ista = Lag+Kmon+Leadtime+1   ! Start index of Y-variable's elements (Initial Forecast Month)
                                ! Tre Lag thang SL so voi thoi diem hien tai, Lay Kmon thang SL tinh tu Lag ve truoc
                                ! Lag=0: Thang ngay truoc thang dung lam du bao
                                ! Leadtime=1: Tinh tu thang tiep theo ngay sau thang dung lam du bao  
   Nspl = Ntime - ista + 1      ! Dung luong mau
   Allocate (X(Kmon*NEigenMax), Y(Ntime))
 !  Return Yobs(Nst,Ntime,NTC_src)
   Do isrc=1, NTC_src
      Do ist=1, Nst
         Do it=1, Ntime
            Y(it)=Yobs(ist,it,isrc)
         Enddo
         it=0
         Do iei=1, NEigenMax
            Do k=1, Kmon        ! Number of Months to be used as predictors
               it=it+1; write(EOFName(it),"(A3,I1,A1,I1)") "EOF",iei,"_",k
            Enddo
         Enddo
         !print "(12A8)", EOFName(:)
         
!         write (Fname,"(A,A,I3.3,A)") trim(Out_Dir)//trim(Yname(isrc)),"_Stn_",ist,".txt"
         write (Fname0,"(A,A,I3.3,A)") trim(Yname(isrc)),"_Stn_",ist
         Fname  = trim(Out_Dir)//trim(Fname0)//".txt"
         Fname1 = trim(Out_Dir)//"Ann_Out_"//trim(Fname0)//".txt"
         Fname2 = trim(Out_Dir)//"Ann_Save_"//trim(Fname0)//".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,I2.2,A,I4.4,A)") "Drought index of "//trim(Yname(isrc))//" at ",Leadtime, &
                                   " months in advance since ", ista-Leadtime,"/",Nyr1, " and Predictors"
         write (3,"(3I6)") Nyr1, Nspl, Mvar
         write (3,"(20I3)") (IDX(i), i=1,Mvar)
         write (3,"(20A15)") trim(Yname(isrc)), (EOFName(k), k=1, Mvar-1)
         Do itime=1, Nspl
            it=0
            Do iei=1, NEigenMax
               Do k=1, Kmon
                  it=it+1
                  X(it) = EOFs(itime+k-1,iei)
                 ! print*, itime+k-1, ista+itime-1
               Enddo
            Enddo
           ! write (3,"(A15,20F15.3)") Sample(ista+itime-1),
           ! TC_Count(ista+itime-1,isrc), (X(it), it=1, Mvar-1)
            write (3,"(20F15.3)") Y(ista+itime-1), (X(it), it=1, Mvar-1)
         Enddo
         Close (3)
      Enddo
   Enddo
   Close (33); Close (44); Close (55)
End Subroutine Create_file_for_Statistical_Models
!----------------------------------------------------------------------------------!
Subroutine Read_Predictand_data (Fname)
  Integer :: iidx, iyr, imon, year, it, it0
  Character (Len=*) :: Fname
  Real, Allocatable :: X(:,:,:),Xlon(:),Xlat(:),Xtime(:)
  Integer :: Nlon, Nlat, Nt

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

     Call rd_ncdf_3d (trim(Fname), Xlon, Xlat, Xtime, Nlon, Nlat, Nt, X, trim(Yname(iidx)))
     if (.not.Allocated (Yobs)) then 
        Nst = Nlon; Ntime = Nt
        Allocate (Yobs(Nst,Ntime,NTC_src))
        Deallocate (Xlon)
        Call rd_ncdf_1d (trim(Fname), Xlon, Nlon, Lon, "lon")
        Deallocate (Xlon)
        Call rd_ncdf_1d (trim(Fname), Xlon, Nlon, Lat, "lat")
        Deallocate (Xlon)
        Call rd_ncdf_1d (trim(Fname), Xlon, Nlon, Reg, "reg")
     endif 
     Yobs(1:Nst,1:Ntime,iidx) = X(1:Nst,1,1:Ntime)
  Enddo
End Subroutine Read_Predictand_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

  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)))
     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_Timeseris_Data !(Fname)
  Integer :: it, iyr, i, j

     it=0
     Do iyr=1,Nyear
        Do imon=1,Nmon
           it=it+1
           write(Sample(it),'(I4.4,I2.2)') YY(iyr), imon
print*,it, iyr,imon, Ntr,Sample(it),ClimIDX(iyr,imon,1:Ntr)
           XX(it,1:Ntr) = ClimIDX(iyr,imon,1:Ntr)
        Enddo
     Enddo
  
     Nth0 = NEigenMax
     print*, it, Ntr, Nth0
     CON = 1.
     Call PTNT(XX,NMAX,MMAX,Ntime,Ntr,CON, Nth0, &
               XMean,Stdev,Corr,K_Eigen,TriRieng,VTORieng,TaiTrong,NhanTo,TT_DAQUAY,NT_DAQUAY)
   ! For other models use:
     EOFs(1:Ntime,1:NEigenMax) = NhanTo(1:Ntime,1:NEigenMax)
   ! For Forecast Models:
     PCs (1:Ntr,1:NEigenMax) = TaiTrong(1:Ntr,1:NEigenMax)
     Eig(1:NEigenMax)        = TriRieng(1:NEigenMax)
     Mean(1:Ntr)             = XMean(1:Ntr)
     Std(1:Ntr)              = Stdev(1:Ntr)
  !Enddo  ! k=1,Nmon
End Subroutine PCA_for_Timeseris_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"