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 :: Nx=160, Ny=80 !Nx=481, Ny=161 
Integer, Parameter :: NMAX=150, MMAX=Nx*Ny
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  :: GPI_Org(:,:,:)   ! Lon,Lat,T
Real, Allocatable  :: GPI(:,:)         ! T, Lon*Lat
!----------------------------------------------------------------------------
Real, Allocatable :: EOFs(:,:), PCs(:,:), Eig(:)
Real, Allocatable :: Mean(:), Std(:), Lon(:), Lat(:), Xtime(:)
Character (Len = 15), Allocatable  :: VarName(:), Sample(:), Yname(:)

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, Fname1, Out_Dir, Out_Dir_Fcs
Integer :: NYr1, NYr2, Nmon, NtrMx, NEigenMax, NTC_src
Integer :: Nyear, Ntime, M_grd
Integer :: Kmon, Lag

Integer :: Nlon, Nlat, Nx1, Ny1, Nxy
Character (Len=100) :: ClimFname, TCFname, PCA_out_Fname

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

  NEigenMax = 5    ! test only
  ClimFname="gpi_new_201308.nc"
  Call Read_GPI(trim(ClimFname)) ; print*, "Call Read_Clim_IDX ... Done!"

  M_grd = Nlon*Nlat
  
  print*, M_grd, Ntime, MMAX
  print*, Lon
  print*, Lat
  pause
!  print*, GPI_Org

  Call Calc_PCA_EOF 
  print*, "PCA Calculation ... Done!"

  Fname="PCAs.nc"; Fname1="EOFs.nc"
  Call Write_PCA_EOF (trim(Fname),trim(Fname1)) 

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

CONTAINS
!----------------------------------------------------------------------------------!
Subroutine Read_GPI (Fname)
   Integer :: ix, iy, iyr, ilon, ilat, it, it0, ilon1, ilon2, ilat1, ilat2
   Character (Len=*) :: Fname
   Real, Allocatable :: Xlon(:), Xlat(:)

   Call rd_ncdf_3d (Fname, lon, lat, xtime, nlon, nlat, ntime, GPI_Org, "gpi")

   Do ix=1, Nlon
      if (lon(ix)>100.and.lon(ix)<101) ilon1=ix
      if (lon(ix)>179.and.lon(ix)<180) ilon2=ix
   Enddo
   Do iy=1, Nlat
      if (lat(iy)>39.and.lat(iy)<40) ilat1=iy
      if (lat(iy)>0.and.lat(iy)<1)   ilat2=iy
   Enddo
   Nx1=ilon2-ilon1+1; Ny1=ilat2-ilat1+1
   Allocate(Xlon(Nx1),Xlat(Ny1))
   ilon=0
   Do ix=ilon1, ilon2 
      ilon=ilon+1; Xlon(ilon)=Lon(ix)
   Enddo
   ilat=0
   Do iy=ilat2, ilat1, -1
      ilat=ilat+1; Xlat(ilat)=Lat(iy)
   Enddo
   Deallocate (Lon, Lat); Allocate (Lon(Nx1),Lat(Ny1))
   Lon(:)=Xlon(:); Lat(:)=Xlat(:)
   Deallocate (XLon, XLat)
   
   Nxy = Nx1*Ny1 ! Nlon*Nlat
   Nlon=Nx1; Nlat=Ny1
   !print*,Nx1, Ny1,ilon1,ilon2,ilat1, ilat2; pause

   Allocate (GPI(Ntime,Nxy))
   it=0
   Do iy=ilat2, ilat1, -1 !1, Nlat
      Do ix=ilon1, ilon2  !1, Nlon
         it=it+1; GPI(1:Ntime,it)=GPI_Org(ix,iy,1:Ntime)
      Enddo
   Enddo
   !print*,GPI; pause
End Subroutine Read_GPI
!-----------------------------------------------------------------------
Subroutine Calc_PCA_EOF !(Fname)
Integer :: iyr, i, j
     CON = 1.
   !----------------- cac bien ------------
     Allocate (Sample(Ntime))
     Do iyr=1, Ntime ! Nyear
        write(Sample(iyr),'(I4)') iyr ! YY(iyr)
        XX(iyr,1:M_grd) = GPI(iyr,1:M_grd)
     Enddo
     Nth0 = NEigenMax
     print*, Ntime, M_grd, Nth0
     Call PTNT(XX,NMAX,MMAX,Ntime,M_grd,CON, Nth0, & 
               XMean,Stdev,Corr,K_Eigen,TriRieng,VTORieng,TaiTrong,NhanTo,TT_DAQUAY,NT_DAQUAY)
     NEigenMax = Nth0
   ! For other models use:
     Allocate (EOFs(1:Ntime,1:NEigenMax))
     EOFs(1:Ntime,1:NEigenMax) = NhanTo(1:Ntime,1:NEigenMax)
   ! For Forecast Models:
     Allocate (PCs (1:M_grd,1:NEigenMax),Eig(1:NEigenMax),Mean(1:M_grd),Std(1:M_grd))
     PCs (1:M_grd,1:NEigenMax) = TaiTrong(1:M_grd,1:NEigenMax)
     Eig(1:NEigenMax)          = TriRieng(1:NEigenMax)
     Mean(1:M_grd)             = XMean(1:M_grd)
     Std(1:M_grd)              = Stdev(1:M_grd)
!!  Enddo  ! k=1,Nmon
End Subroutine Calc_PCA_EOF 
!-----------------------------------------------------------------------
Subroutine Write_PCA_EOF (Fnm_PCA, Fnm_EOF)
   Character (Len=*) :: Fnm_PCA, Fnm_EOF
   Real, Allocatable :: X(:,:,:), Ymean(:,:), Ystd(:,:)
   Integer :: ilon,ilat,ik,it

   Allocate (X(1:Nlon,1:Nlat,1:NEigenMax))! ,Y(1:Nlon,1:Nlat))
   it=0
   Do ilat=1, Nlat
      Do ilon=1, Nlon
         it=it+1; X(ilon,ilat,1:NEigenMax)=PCs(it,1:NEigenMax)
      Enddo
   Enddo
   Call wr_ncdf_3d (Fnm_PCA,lon,lat,Eig,nlon,nlat,NEigenMax,X,"PCA","hours since 2010-08-01 ", &
                    "PCA_GPI_6hrs", "-", "lon", "lat","Eig", Miss)
   Deallocate (X)
   Allocate (Ymean(1:Nlon,1:Nlat), Ystd(1:Nlon,1:Nlat))
   it=0
   Do ilat=1, Nlat
      Do ilon=1, Nlon
         it=it+1; Ymean(ilon,ilat)=Mean(it); Ystd(ilon,ilat)=Std(it)
      Enddo
   Enddo
   Call wr_ncdf_2d_apd (Fnm_PCA, Ymean, "Mean",nlon,nlat,"GPI mean", "-", "lon", "lat", Miss)
   Call wr_ncdf_2d_apd (Fnm_PCA, Ystd,  "std" ,nlon,nlat,"GPI_Std",  "-", "lon", "lat", Miss)
   Deallocate (Ymean,Ystd)
   Call wr_ncdf_2d (Fnm_EOF,Xtime,Eig,Ntime,NEigenMax, EOFs, "EOF", "GPI_EOF", "-", "time", "Eig", Miss)    
End Subroutine Write_PCA_EOF
!-----------------------------------------------------------------------
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"