Program FA_ENSO_Monthly
Use NetCDF
Integer, Parameter :: NMAX=1000, MMAX=20
Real, Parameter    ::    Miss=-99.

Real               :: XX(Nmax, Mmax)
Integer            :: YY(Nmax)

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

Character (Len = 15), Allocatable  :: VarName(:), Sample(:), Yname(:)
Real, Allocatable    :: ClimIDX(:,:,:)
Real, Allocatable    :: YTCs(:,:,:)
Integer, Allocatable :: mon_x(:), yrs_x(:)
Character (Len=3), Allocatable :: MonName(:)

Real, Allocatable    :: EOFs(:,:,:) 

Integer :: Nth0     ! So gia tri rieng duoc giu lai
!----------------------------------------------------------------------------
Integer :: NTram, IO, yr, Ntr, iyr, ist, imon, K_Eigen
Character (Len = 255) :: Fname, Out_Dir, Out_Dir_Fcs
Character (Len = 100) :: ClimFname, TCFname, PCA_out_Fname

Integer :: NYr1, NYr2, Nmon, NtrMx, NEigenMax, NTC_src
Integer :: Nyear, Ntime
Integer :: Kmon, Lag

Integer :: mon_y
Integer :: Nyrs, yrs_y

Character (Len=4) :: Var
Real    :: Con

Namelist /Lag_Kmon/ Lag, Kmon
Namelist /Data_info_Fcs/ 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_Fcs); 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)) 
    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 

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

  print*, "Reading Input data... Done!"

  Call Calc_EOFs
  Call Create_file_for_Statistical_Models (Lag, Kmon)
  Call Make_Namelist_for_TC_Traing_models

  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_Fcs.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 =    .F.       ! 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 Calc_EOFs
  Real, Allocatable :: PCs(:,:,:), Mean(:,:), Std(:,:), Eig(:,:), Eig1(:), Eig2(:), Eig3(:), Nspl(:) 
  Real, Allocatable :: XMon(:), XIDX(:), XEig(:)
  Character (Len=20) :: var_nm
  Integer :: nlon,nlat,nt, Nsample  ! Tuong ung la Nmon, Ntr, NEigMx
  Real :: Xorg(1:Nyear,1:NtrMx), Xpre(NtrMx), Fpre(NEigenMax), A(NtrMx,NEigenMax), Tmp 

  ! Read all PCs, Mean, Std and Eigen values from previous calculations:
  Allocate (Eig(Nmon,NEigenMax))
  Call rd_ncdf_3d ("PCs_and_Eig_Values.nc",XMon,XIDX,XEig,nlon,nlat,nt, PCs, "PC")
  Deallocate (XMon,XIDX,XEig)
  Call rd_ncdf_2d ("PCs_and_Eig_Values.nc",XMon,XIDX,nlon,nlat, Mean, "Mean")
  Deallocate (XMon,XIDX)
  Call rd_ncdf_2d ("PCs_and_Eig_Values.nc",XMon,XIDX,nlon,nlat, Std, "Std")
  Deallocate (XMon,XIDX)
  Call rd_ncdf_1d ("PCs_and_Eig_Values.nc",XMon,nlon, Eig1, "Eig1")
  Deallocate (XMon)
  Call rd_ncdf_1d ("PCs_and_Eig_Values.nc",XMon,nlon, Eig2, "Eig2")
  Deallocate (XMon)
  Call rd_ncdf_1d ("PCs_and_Eig_Values.nc",XMon,nlon, Eig3, "Eig3")
  Deallocate (XMon)
  Call rd_ncdf_1d ("PCs_and_Eig_Values.nc",XMon,nlon, Nspl, "Nspl")

  Eig(:,1)=Eig1(:); Eig(:,2)=Eig2(:); Eig(:,3)=Eig3(:)
  Deallocate (XMon,Eig1,Eig2,Eig3)
  ! Calculating the EOFs (F): X(N,M) x A(M,K) = F(N,K)
  Do imon=1, Nmon
     tmp=Sqrt(Nspl(imon))
     A(1:Ntr,1:NEigenMax)=PCs(imon,1:Ntr,1:NEigenMax)    ! Cac PCs
     Do iyr=1,Nyear
        write(Sample(iyr),'(I4)') YY(iyr)
        Xorg(iyr,1:Ntr) = ClimIDX(iyr,imon,1:Ntr)
        Do j=1, Ntr; Xpre(j) = (Xorg(iyr,j)-Mean(imon,j))/Std(imon,j)/tmp; Enddo
        Fpre=Matmul(Xpre,A)
        Do j=1, NEigenMax; Fpre(j)=Fpre(j)/Eig(imon,j); Enddo
        EOFs(iyr,1:NEigenMax,imon)=Fpre(1:NEigenMax)
     Enddo
  Enddo

End Subroutine Calc_EOFs
!-----------------------------------------------------------------------
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

   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 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_Fcs)//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
      print "(I4,12F8.3)", iyr, YTCs(it0,:,iidx)
   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
      print "(I4,12F8.3)", iyr, ClimIDX(it0,:,iidx)
   Enddo
Enddo
End Subroutine Read_Clim_IDX
!-----------------------------------------------------------------------
Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"

END  ! Main Prog
!-----------------------------------------------------------------------