Use netCDF
      ! Chuong trinh tinh he so tuong quan giua So luong TC tren VES va cac truong KQ-DD cua CFS_Rfc
      ! Chuoi thoi gian tu 1982-2009
      ! Lead time lay tu 1-6 thang
      ! TC co the la monthly hoac seasonally (tc, tc3m,...). Co the thay TC bang NTC
      ! Cac truong co the la monthly hoac seasonally. Neu la seasonally thi lay trung binh truot
      ! Mien tinh la khu vuc TBD-ADD duoc xac dinh boi Lon1, Lon2, Lat1, Lat2 
      Integer, Parameter :: Nmon=12, Yr1=1982, Yr2=2009, Nyrs=(Yr2-Yr1)+1
      Integer, Parameter :: Ntime=Nyrs*Nmon, LT=6
      Real, Parameter :: Miss=-99.0, missing_value = 9.999e+20
      Real, Parameter :: Dlon=1.0, Dlat=1.0 !, Lon1=60.0, Lon2=290.0, Lat1=-20.0, lat2=60.0
      Character (Len=11), Parameter :: Period=".1982-2009."
      Character (Len=15), Parameter :: Subdir="ensmean_all/"

      Integer :: Ntag_mon, Nst,Nlon,Nlat, ilon,ilat,ist, ilt, NLT, jlt, iv, lonlat

      Real, Allocatable :: sst(:,:,:,:,:)    ! (Nlon,Nlat,Nyrs,Nmon,LT)
      Real, Allocatable :: Xall(:,:,:,:,:,:) ! (Nlon,Nlat,Nyrs,Nmon,LT,Nvar)

      Real, Allocatable :: TC(:,:)        ! (Nyrs,Nmon)
      Real, Allocatable :: xlon(:),xlat(:),xyrs(:),xmon(:)  ! Nlon,Nlat,Nyrs,Nmon
      Character (Len=100) :: Fname, Dir
      Character (Len=10)  :: vname
      Character (Len=1)   :: Lst
      !-------------------------------------------------------------------------------------------
      Real, Allocatable :: X1D(:), X2D(:,:), X3D(:,:,:)  ! (Nlon,Nlat,Nt)
      !-------------------------------------------------------------------------------------------
      Integer, Parameter :: Nvar=8
      Integer :: lev, ivar
      ! AO fields:
      Character (Len=5) :: Vars(Nvar)= &
                (/"sst","prmsl","vor","vws","vvel","rh","h85","h500"/)

      !         (/"sst","vor","vws","vvel","rh","h85"/) !,"h500"/)
      !         (/"sst","prmsl","vor","vws","vvel","rh","h85"/) !,"h500"/)
      !                                     850,   700, 850
      !Character (Len=10), Parameter :: vartc="ntc3m"
      Character (Len=10) :: vartc
      Integer, Parameter :: Nvartc=4
      Character (Len=10), Parameter :: vartcn(Nvartc)= &
                           (/"tc3m","tc6m","ntc3m","ntc6m"/)

      Character (Len=20), Allocatable :: Varn_out(:,:,:)
      Integer :: NSpl(Nvar)   ! Dung luong mau
      !Integer, Allocatable :: Targ_mon(:), Mvar(:,:), Mvung(:,:,:), LeadTime(:)
      Real, Allocatable :: Vung(:,:), Xfc(:,:,:,:)  ! (LT, TarMon, Nyrs, Nxy): So lieu cho mo hinh du bao
      Integer :: ii, jj, kk, Nxy, Nxmx, Nymx
      Integer, Allocatable :: Nxpc(:), Nypc(:) 
      Real, Allocatable    :: lonpc(:,:), latpc(:,:) ! Dung cho PCA
      !-------------------------------------------------------------------------------------------
      ! For PCA:
      Real*8, Allocatable :: xdat(:,:), ydat(:,:), xtb(:), sx(:)
      Real*8, Allocatable :: EOFs(:,:), PCs(:,:)

      Integer :: Nt, Nx   ! Kich thuoc thuc cua mang X
      Integer :: Mcount, Ncount, switch, Standard, PC_Stard, Nmode, K
      Real*8, Allocatable   :: Eigenvalues(:), ECumulate(:)

      Real*8 :: Factor

      ! Write out for calculating EOFs:
      Real, Allocatable :: PCall(:,:,:,:)              ! (LT,Nmon,Nx,Nmode)
      Real, Allocatable :: EOFall(:,:,:,:)             ! (LT,Nmon,Nx,Nmode)
      Real, Allocatable :: xtball(:,:,:), sxall(:,:,:) ! (LT,Nmon,Nx)
      Real, Allocatable :: Eigall(:,:,:), Nsple(:,:,:)  ! (LT,Nmon,Nx)
      Real, Allocatable :: EVari (:,:,:), ECum (:,:,:)  ! (LT,Nmon,Nx)

      Character (Len=200) :: pca_out, eof_out, Out_Dir
      
      !-------------------------------------------------------------------------------------------
      Open (1,file="Extract_data_New.txt", status="old")
      Allocate (Vung(Nvar,4))  ! 4=lon1,lon2,lat1,lat2
      read (1,*)
      Do ivar=1, Nvar
         read (1,*) ii, (Vung(ivar,jj),jj=1,4)
      Enddo
      Close (1)

      !Fname="/home/tanpv/Data/TCs/VES_Monthly_TC_NTC_Count.nc"
      !-------------------------------------------------------------------------------------------
      !vartc="tc3m"
      !-------------------------------------------------------------------------------------------
      ! Reading TC data:
      !Call Read_TC (Fname,trim(vartc))
      
      ! Reading AO fields:
      Dir="/home/tanpv/Progs/TCs/Data/CFS_Rfc/"//trim(Subdir) !ensmean_all/"
      !---------------------------------------------------------------------------------
      ! "sst","prmsl","vor","vws","vvel","rh","h85","h500"
      !   1      2      3     4      5     6    7      8
      Do ivar=1,Nvar
         vname=trim(Vars(ivar))
         print*,"Calculating for ",vname
         !if (ivar <= 3) then      ! Reading 3D-arrays (lon,lat,time)
         !   Call Read_var (vname)
         !else if (ivar<=8) then   ! Reading 4D-arrays (lon,lat,lev,time)
            if (trim(Vars(ivar))=="sst".OR.trim(Vars(ivar))=="prmsl".OR. &
                trim(Vars(ivar))=="vor".OR.trim(Vars(ivar))=="vws") then
                Call Read_var (vname)
            else if (trim(Vars(ivar))=="h85") then
               lev=850
               Call Read_var_4D ("hgt",lev)
            else if (trim(Vars(ivar))=="h500") then
               lev=500
               Call Read_var_4D ("hgt",lev)
            else if (trim(Vars(ivar))=="vvel") then
               lev=850
               Call Read_var_4D (vname,lev)
            else if (trim(Vars(ivar))=="rh") then
               lev=700
               Call Read_var_4D (vname,lev)
            endif
         !endif
         if (ivar==1) Allocate(Xall(Nlon,Nlat,Nyrs,Nmon,LT,Nvar))
         Xall(:,:,:,:,:,ivar)=sst(:,:,:,:,:)
         print*,trim(Vars(ivar)), Nlon,Nlat,Nyrs,Nmon,LT,Nvar
         !print*,Xall(130,110,3:5,7,3,ivar)
         !Xall(:,:,:,:,:,ivar)=sst(:,:,:,:,:)

         !Do ilt=1,LT
         !   Do iyr=1,Nyrs
         !      print"(I5,12F10.3)",iyr+Yr1-1,sst(20,20,iyr,1:Nmon,ilt)
         !   Enddo
         !   pause
         !Enddo
      Enddo
      print*,"Reading all data ... Done !"
      Call Extract_Var_1

      !---------------------------------------------------------------------------------
      ! PCA:
      !Open (1,file="Namelist_PCA_New.txt",Status="old")
      !Read (1, Model_Info)
      !Close (1)
      !---------------------------------------------------------------------------------
      Nt = Nyrs-1; Nx = Nxy
      switch     =  1            ! =1: Use X(nx,nt) instead of X(nt,nx) ; =0: Origin X(nt,nx)
      standard   =  1           ! =1: Su dung ma tran Correlation;       =0: Su dung ma tran Covariance
      PC_Stard   =  0         ! =1: PCs chuan hoa (PCs=EVtors)         =0: EOFs chuan hoa
      Nmode      =  6         ! So mode duoc giu lai (Nmode <= Min{nt,nx})

      !=======================================      
      ! Test with Timeseries=1983-2002 (Nt=20)
      Nt = 20     ! De lai 7 nam (2003-2009)
      !=======================================      

      if (switch==1) then
         Mcount=Nt; Ncount=Nx
      else if (switch==0) then
         Mcount=Nx; Ncount=Nt
      endif
      Factor = Sqrt(Real(Nt))

      Allocate (xdat(Nt,Nx), ydat(Nt,Nx))
      Allocate (Xtb(Nx), Sx(Nx))   ! TB, DLC, ma tran Covar hoac Corre
      Allocate (Eigenvalues(Mcount))          ! Tri rieng va Vecto rieng (1D)
      Allocate (EOFs(Nt,Nmode), PCs(Nx,Nmode), ECumulate(Nmode))

      Allocate (PCall(LT,Nmon,Nx,Nmode), EOFall(LT,Nmon,Nt,Nmode))
      Allocate (xtball(LT,Nmon,Nx), sxall(LT,Nmon,Nx))
      Allocate (Eigall(LT,Nmon,Nx), Nsple(LT,Nmon,Nx))
      Allocate (EVari(LT,Nmon,Nx), ECum(LT,Nmon,Nx))
      PCall=Miss; EOFall=Miss; xtball=Miss; sxall=Miss; Eigall=Miss; Nsple=Miss
      EVari=Miss; ECum=Miss

      !---------------------------------------------------------------------------------
      Do ilt=1,LT
         Do imon=1,Nmon
            Do it=1,Nt
               xdat(it,1:Nx)=Xfc(ilt,imon,it,1:Nx)
            Enddo
            Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, &
                           Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode)
            Do ivar=1,Nvar
               write (pca_out,"(A,I1,A,I2.2,A)") "Dat_out/pca_tc3m_"//trim(Vars(ivar))//"_LT_",ilt,"_Mon_",imon,".nc"
               write (eof_out,"(A,I1,A,I2.2,A)") "Dat_out/eof_tc3m_"//trim(Vars(ivar))//"_LT_",ilt,"_Mon_",imon,".nc"
            !   Call Write_EOF_PCA_for_Plot (pca_out,eof_out,ivar)
            Enddo
            PCall (ilt,imon,1:Nx,1:Nmode) = PCs (1:Nx,1:Nmode)
            EOFall(ilt,imon,1:Nt,1:Nmode) = EOFs(1:Nt,1:Nmode)
            xtball(ilt,imon,1:Nx)         = Xtb (1:Nx)
            sxall (ilt,imon,1:Nx)         = Sx  (1:Nx)
            Eigall(ilt,imon,1:Mcount)     = Eigenvalues(1:Mcount)
            Nsple (ilt,imon,1:Nx)         = Nt
            ! Test:
            !print"(5F12.3)",Matmul(Transpose(PCs),PCs) !==> Lamda
            !print*
            !print"(5F12.3)",Matmul(Transpose(EOFs),EOFs) !==> 1
            !pause
         Enddo
      Enddo 
      ! Write out for EOF calculation:
      !pca_out="Dat_out/pca_tc3m_all.nc" 
      !eof_out="Dat_out/eof_tc3m_all.nc" 
      !Call Write_EOF_PCA_for_Fcs (pca_out,eof_out)

      Out_Dir="/home/tanpv/Progs/TCs/FCS_Aug/Out_train_eof_cfs_new/"
      Do jtc=1,Nvartc
         Fname="/home/tanpv/Data/TCs/VES_Monthly_TC_NTC_Count_jp_new.nc"
         vartc=vartcn(jtc)
         Call Read_TC (Fname,trim(vartc))

         ilt=3
         !Out_Dir="/home/tanpv/Progs/TCs/FCS/Out1/"
         Call Create_file_for_Statistical_Models(lt,Out_Dir)

         Deallocate (TC)
      Enddo
      !---------------------------------------------------------------------------------
      ! Sl Doc lap:
      Deallocate (EOFall)
      Nt = Nyrs-1
      Call Calc_EOFs_for_Fcs_DL

      Out_Dir="/home/tanpv/Progs/TCs/FCS_Aug/Out_fcs_eof_cfs_new/"
      Do jtc=1,Nvartc
         Fname="/home/tanpv/Data/TCs/VES_Monthly_TC_NTC_Count_jp_new.nc"
         vartc=vartcn(jtc)
         Call Read_TC (Fname,trim(vartc))

         !Nt = Nyrs-1
         !Call Calc_EOFs_for_Fcs_DL
         ilt=3
         !Out_Dir="/home/tanpv/Progs/TCs/FCS/Out_Fcs1/"
         Call Create_file_for_Fcs_DLap (ilt,Out_Dir)
         Deallocate (TC)
      Enddo

      !---------------------------------------------------------------------------------
      CONTAINS
      !-------------------------------------------------------------------------------------------
      Subroutine Extract_Var_1
        Integer :: it,ilt,imon,ilon,ilat,ix,iy,ivar,iv1,iv2,jlt, Nspl_ivar
        Real :: xlon1,xlon2,xlat1,xlat2
        Integer :: NSpl1(Nvar), Nxy1

        Call Get_Domain_Size_and_Lon_Lat  ! Return: Nxy, lonpc, latpc
        print*,"Begining Extract_Var_1,LT,Nmon,Nyrs,Nxy,Nlon,Nlat,Nyrs=",LT,Nmon,Nyrs,Nxy,Nlon,Nlat,Nyrs
        Allocate (X3D(Nlon,Nlat,Nyrs))
        Allocate (Xfc(LT,Nmon,Nyrs-1,Nxy))

        Do ilt=1,LT
           Do imon=1,Nmon 
              iv1=0
              Do ivar=1,Nvar
                 !print*,"ilt,imon,ivar=",ilt,imon,ivar; !pause
                 Nspl_ivar=NSpl(ivar); !iv2=iv1+Nspl_ivar-1
                 Allocate (X2D(Nyrs,Nspl_ivar))
                 X3D(1:Nlon,1:Nlat,1:Nyrs)=Xall(1:Nlon,1:Nlat,1:Nyrs,imon,ilt,ivar)
                 xlon1=Vung(ivar,1); xlon2=Vung(ivar,2)
                 xlat1=Vung(ivar,3); xlat2=Vung(ivar,4)
                 Call Subdomain_1 (xlon1,xlon2,xlat1,xlat2)  
                 !print*,"iv1,iv2,Nspl_ivar=",iv1,iv2, Nspl_ivar, iv2-iv1+1
                 Do it=1,Nspl_ivar
                    iv1=iv1+1
                    !print*,iv1,it,Nspl_ivar
                    Do iyr=1,Nyrs-1
                       !Xfc(ilt,imon,1:Nyrs,iv1)=X2D(1:Nyrs,it)
                       Xfc(ilt,imon,iyr,iv1)=X2D(iyr+1,it)
                       !if (Xfc(ilt,imon,iyr,iv1)==Miss) print*,ivar,iv1,lonpc(iv1),latpc(iv1)
                    Enddo
                 Enddo
                 !print*,"iv1,Nspl_ivar=",iv1,Nspl_ivar
                 !iv1=iv2+1
                 Deallocate (X2D)
              Enddo
           Enddo
        Enddo
      End Subroutine Extract_Var_1
      !-------------------------------------------------------------------------------------------
      Subroutine Get_Domain_Size_and_Lon_Lat !(ilt,itarg)
         Integer :: ilo, ila, iyr, ilon,ilat,ivar, Nx,Ny !, Nxy
         Integer :: ilt,itarg    ! Leadtime, Target month
         Real :: xlon1,xlon2,xlat1,xlat2

         Allocate (Nxpc(Nvar), Nypc(Nvar))
         Nxy=0; Nxmx=0; Nymx=0
         Do ivar=1,Nvar
            xlon1=Vung(ivar,1); xlon2=Vung(ivar,2); 
            xlat1=Vung(ivar,3); xlat2=Vung(ivar,4);
            Nx=Int((xlon2-xlon1)/Dlon)+1; Ny=Int((xlat2-xlat1)/Dlat)+1
            Nxpc(ivar)=Nx; Nypc(ivar)=Ny
            if (Nxmx<Nx) Nxmx=Nx; if (Nymx<Ny) Nymx=Ny     ! Max size of Nlon, Nlat
            NSpl(ivar)=Nx*Ny  ! Dung luong mau cua ivar
            Nxy=Nxy+Nx*Ny     ! Dung luong mau cua Nvar
         Enddo
         !print*,"Samples: ",NSpl, Nxy; ! pause
         ! 
         Allocate (lonpc(Nvar,Nxmx),latpc(Nvar,Nymx))

         Do ivar=1,Nvar
            xlon1=Vung(ivar,1); xlon2=Vung(ivar,2);
            xlat1=Vung(ivar,3); xlat2=Vung(ivar,4);
            Nx=Int((xlon2-xlon1)/Dlon)+1; Ny=Int((xlat2-xlat1)/Dlat)+1
            ix=Int(xlon1/Dlon)+1; iy=Int((xlat1+90)/Dlat)+1
            ila=0
            Do ilat=iy, iy+Ny-1
               ila=ila+1; latpc(ivar,ila)=xlat(ilat)
               ilo=0
               Do ilon=ix, ix+Nx-1
                  ilo=ilo+1; lonpc(ivar,ilo)=xlon(ilon)
               Enddo
            Enddo
         Enddo
      End Subroutine Get_Domain_Size_and_Lon_Lat
      !-------------------------------------------------------------------------------------------
      Subroutine Subdomain_1 (xlon1,xlon2,xlat1,xlat2)
         Real :: xmean,xlon1,xlon2,xlat1,xlat2, tmp
         Integer :: iyr, ilon,ilat,ic, Nx,Ny
         Real, Allocatable :: X2D1(:,:) !(Nlon,Nlat)

         Nx=Int((xlon2-xlon1)/Dlon)+1; Ny=Int((xlat2-xlat1)/Dlat)+1
         ix=Int(xlon1/Dlon)+1; iy=Int((xlat1+90)/Dlat)+1

         ! X2D(Nyrs,Nxy),X3D(Nlon,Nlat,Nyrs))
         Allocate (X2D1(Nlon,Nlat))
         Do iyr=1, Nyrs
            X2D1(1:Nlon,1:Nlat)=X3D(1:Nlon,1:Nlat,iyr)
            ic=0
            Do ilat=iy, iy+Ny-1 ! ilon=ix, ix+Nx
               Do ilon=ix, ix+Nx-1 !ilat=iy, iy+Ny
                  ic=ic+1; X2D(iyr,ic)=X2D1(ilon,ilat)
               Enddo
            Enddo
            !print*,iyr,ic
         Enddo
         Deallocate (X2D1)
      End Subroutine Subdomain_1

      !-------------------------------------------------------------------------------------------
      Subroutine Read_TC (Fname,vname)
        ! vname: tc, tc3m,...
        Character (Len=*) :: Fname, vname
        Real, Allocatable :: xyear(:),xmon(:), xtc(:,:)
        Integer :: ilt, Ny, Nm, Nt, it, ix, iy, imon, iyr

        Call rd_ncdf_2d (Fname,xyear,xmon,Ny,Nm,xtc,trim(vname))
        Allocate (TC(Nyrs,Nmon))
        iy=0
        Do iyr=Yr1,Yr2   ! Chi lay chuoi thoi gian tu Yr1-Yr2
           Do it=1,Ny
              if (int(xyear(it))==iyr) then
                 iy=iy+1; TC(iy,1:Nmon)=xtc(it,1:Nmon)
              endif
           Enddo
           !print"(I5,12F5.1)",iy,(TC(iy,1:Nmon))
        Enddo
        ! Luu vao mang Xfc:
        !Do ilt=1, NLT
        !   Do it=1,Ntag_mon
        !      imon=Targ_mon(it)
        !      Xfc(ilt,it,1:Nyrs,1)=TC(1:Nyrs,imon)  ! Bien phu thuoc
        !   Enddo
        !Enddo
      End Subroutine Read_TC
      !-------------------------------------------------------------------------------------------
      Subroutine Read_var_4D_mlv (vname)
        Character (Len=*) :: vname
        Integer :: Nt, Nx, Ny, Nlev, it, ix, iy, imon, iyr, ilt, lev,ilev
        Integer :: ilv1, ilv2, ilv3
        Real, Allocatable :: xtmp(:,:,:), X(:,:,:,:), lon(:),lat(:),xt(:)
        Real, Allocatable :: lonx1(:),latx1(:),lev1(:),xt1(:),XX(:,:,:,:)
        Real, Allocatable :: X1(:,:,:,:),X2(:,:,:,:),X3(:,:,:,:)
        Real, Allocatable :: xave(:)
        Character (Len=250) :: Fname

        ilv1=Level(700); ilv2=Level(600); ilv3=Level(300)
        Do ilt=1,LT
           write (Lst,"(I1)") ilt
           Fname=trim(Dir)//Lst//"/"//trim(vname)//Period//"lead"//Lst//".nc"
           Call rd_ncdf_4d (trim(Fname),lonx1,latx1,lev1,xt1,nx,ny,nlev,nt,X,trim(vname))
           if (ilt==1) then
              Allocate (X1(nx,ny,Ntime,LT),lon(nx),lat(ny),xt(nt))
              Allocate (X2(nx,ny,Ntime,LT),X3(nx,ny,Ntime,LT))
              Allocate (XX(nx,ny,Ntime,LT))
              lon(:)=lonx1(:); lat(:)=latx1(:); xt(:)=xt1(:)
           endif
           X1(1:Nx,1:Ny,1:Nt,ilt)=X(1:Nx,1:Ny,ilv1,1:Nt)  ! Chi lay muc ilev
           X2(1:Nx,1:Ny,1:Nt,ilt)=X(1:Nx,1:Ny,ilv2,1:Nt)
           X3(1:Nx,1:Ny,1:Nt,ilt)=X(1:Nx,1:Ny,ilv3,1:Nt)
           Deallocate (lonx1,latx1,lev1,xt1,X)
        Enddo
        XX=(X1+X2+X3)/3.0      ! Trung binh 3 truong
        Deallocate (X1, X2, X3)
        Nlon=Nx; Nlat=Ny

        !Nlon=Int((Lon2-Lon1)/Dlon)+1; Nlat=Int((lat2-Lat1)/Dlat)+1  ! Mien TBD-ADD
        !ix=Int(Lon1/Dlon)+1; iy=Int((Lat1+90)/Dlat)+1

        if (.not.Allocated (xlon)) Allocate (xlon(Nlon))
        if (.not.Allocated (xlat)) Allocate (xlat(Nlat))
        Allocate (xtmp(Nlon,Nlat,Ntime))      ! Lat, lon for plotting
        if (.not.Allocated(sst)) Allocate (sst(Nlon,Nlat,Nyrs,Nmon,LT))
        Allocate (xave(Ntime))

        !Allocate (xlon(Nlon),xlat(Nlat),xtmp(Nlon,Nlat,Ntime))      ! Lat, lon for plotting
        !Allocate (sst(Nlon,Nlat,Nyrs,Nmon,LT),xave(Ntime))

        xlon(1:Nlon) = lon(1:Nlon); xlat(1:Nlat) = lat(1:Nlat)
        Do ilt=1,LT
           xtmp(:,:,:)=Miss ! missing_value
           xtmp(1:Nlon,1:Nlat,ilt+1:Ntime) = XX(1:Nlon,1:Nlat,1:Nt,ilt)
           Where (xtmp==missing_value) xtmp=Miss     ! Convert to my missing values
           it=0
           Do iyr=1,Nyrs       ! Resharp to 5D-array (lon,lat,year,mon,leadtime)
              Do imon=1,Nmon
                 it=it+1
                 sst(1:Nlon,1:Nlat,iyr,imon,ilt)=xtmp(1:Nlon,1:Nlat,it)
              Enddo
           Enddo
        Enddo
        !Where (sst==missing_value) sst=Miss
        Deallocate (xtmp,lon,lat,xt,XX,xave)
      End Subroutine Read_var_4D_mlv
      !-------------------------------------------------------------------------------------------
      Subroutine Read_var_4D (vname,lev)
        Character (Len=*) :: vname
        Integer :: Nt, Nx, Ny, Nlev, it, ix, iy, imon, iyr, ilt, lev,ilev
        Real, Allocatable :: xtmp(:,:,:), X(:,:,:,:), lon(:),lat(:),xt(:)
        Real, Allocatable :: XX(:,:,:,:), lonx1(:),latx1(:),lev1(:),xt1(:)
        Real, Allocatable :: xave(:)
        Character (Len=250) :: Fname

        ilev=Level(lev) 
        Do ilt=1,LT
           write (Lst,"(I1)") ilt
           Fname=trim(Dir)//Lst//"/"//trim(vname)//Period//"lead"//Lst//".nc"
           Call rd_ncdf_4d (trim(Fname),lonx1,latx1,lev1,xt1,nx,ny,nlev,nt,X,trim(vname))

           if (ilt==1) then
              Allocate (XX(nx,ny,Ntime,LT),lon(nx),lat(ny),xt(nt))
              lon(:)=lonx1(:); lat(:)=latx1(:); xt(:)=xt1(:)
           endif
           XX(1:Nx,1:Ny,1:Nt,ilt)=X(1:Nx,1:Ny,ilev,1:Nt)  ! Chi lay muc ilev
           Deallocate (lonx1,latx1,lev1,xt1,X)
        Enddo
        Nlon=Nx; Nlat=Ny;

        !Nlon=Int((Lon2-Lon1)/Dlon)+1; Nlat=Int((lat2-Lat1)/Dlat)+1  ! Mien TBD-ADD
        !ix=Int(Lon1/Dlon)+1; iy=Int((Lat1+90)/Dlat)+1
        if (.not.Allocated (xlon)) Allocate (xlon(Nlon))
        if (.not.Allocated (xlat)) Allocate (xlat(Nlat))
        Allocate (xtmp(Nlon,Nlat,Ntime))      ! Lat, lon for plotting
        if (.not.Allocated(sst)) Allocate (sst(Nlon,Nlat,Nyrs,Nmon,LT))
        Allocate (xave(Ntime))

        !xlon(1:Nlon) = lon(ix:ix+Nlon); xlat(1:Nlat) = lat(iy:iy+Nlat)
        Do ilt=1,LT
           xtmp(:,:,:)=Miss ! missing_value
           xtmp(1:Nlon,1:Nlat,ilt+1:Ntime) = XX(1:Nlon,1:Nlat,1:Nt,ilt)
           Where (xtmp==missing_value) xtmp=Miss     ! Convert to my missing values
           it=0
           Do iyr=1,Nyrs       ! Resharp to 5D-array (lon,lat,year,mon,leadtime)
              Do imon=1,Nmon
                 it=it+1
                 sst(1:Nlon,1:Nlat,iyr,imon,ilt)=xtmp(1:Nlon,1:Nlat,it)
              Enddo
           Enddo
        Enddo
        !Where (sst==missing_value) sst=Miss
        if (trim(vname)=="vvel") then
           Where (sst/=Miss) sst=sst*100.0      ! x 10^2
        else if (trim(vname)=="hgt") then
           Where (sst/=Miss) sst=sst/10.0       ! / 10 m-->dam
        endif

        Deallocate (xtmp,lon,lat,xt,XX,xave)
      End Subroutine Read_var_4D
      !-------------------------------------------------------------------------------------------
      Integer Function Level (lev)
        Integer :: lev,ilev
        if (lev==1000) then
           ilev=1
        elseif (lev==925) then
           ilev=2
        elseif (lev==850) then
           ilev=3
        elseif (lev==700) then
           ilev=4
        elseif (lev==600) then
           ilev=5
        elseif (lev==500) then
           ilev=6
        elseif (lev==300) then
           ilev=7
        elseif (lev==200) then
           ilev=8
        elseif (lev==100) then
           ilev=9
        endif
        Level=ilev
      End Function Level
      !-------------------------------------------------------------------------------------------
      Subroutine Read_var (vname)
        Character (Len=*) :: vname
        Integer :: Nt, Nx, Ny, it, ix, iy, imon, iyr, ilt
        Real, Allocatable :: xtmp(:,:,:), X(:,:,:), lon(:),lat(:),xt(:)
        Real, Allocatable :: XX(:,:,:,:), lonx1(:),latx1(:),xt1(:)
        Character (Len=250) :: Fname
        Real, Allocatable :: xave(:)
        Real :: tmp

        Do ilt=1,LT
           write (Lst,"(I1)") ilt
           Fname=trim(Dir)//Lst//"/"//trim(vname)//Period//"lead"//Lst//".nc"
           Call rd_ncdf_3d (trim(Fname),lonx1,latx1,xt1,nx,ny,nt, X, trim(vname))
           if (ilt==1) then 
              Allocate (XX(nx,ny,Ntime,LT),lon(nx),lat(ny),xt(nt))
              lon(:)=lonx1(:); lat(:)=latx1(:); xt(:)=xt1(:)
           endif
           XX(1:Nx,1:Ny,1:Nt,ilt)=X(1:Nx,1:Ny,1:Nt)  ! Data from all leadtimes
           Deallocate (lonx1,latx1,xt1,X)
        Enddo
        
        Nlon=Nx; Nlat=Ny;
        if (.not.Allocated (xlon)) Allocate (xlon(Nlon))
        if (.not.Allocated (xlat)) Allocate (xlat(Nlat))
        Allocate (xtmp(Nlon,Nlat,Ntime))      ! Lat, lon for plotting
        if (.not.Allocated(sst)) Allocate (sst(Nlon,Nlat,Nyrs,Nmon,LT))
        Allocate (xave(Ntime))

        xlon(1:Nlon) = lon(1:Nlon); xlat(1:Nlat) = lat(1:Nlat)

        Do ilt=1,LT
           xtmp(:,:,:)=Miss ! missing_value
           xtmp(1:Nlon,1:Nlat,ilt+1:Ntime) = XX(1:Nlon,1:Nlat,1:Nt,ilt)

           Where (xtmp==missing_value) xtmp=Miss
           !
           it=0
           Do iyr=1,Nyrs       ! Resharp to 5D-array (lon,lat,year,mon,leadtime)
              Do imon=1,Nmon
                 it=it+1
                 sst(1:Nlon,1:Nlat,iyr,imon,ilt)=xtmp(1:Nlon,1:Nlat,it)
              Enddo
           Enddo
        Enddo

        if (trim(vname)=="sst") then
           Where (sst/=Miss) sst=sst-273.15 ! Convert to C degree
        else if (trim(vname)=="prmsl") then
           Where (sst/=Miss) sst=sst/100.0
        else if (trim(vname)=="vor") then
           Where (sst/=Miss) sst=sst*1000000.0  ! x 10^6
        endif

        print*,nx,ny,nt,Nlon,Nlat,ntime
        print*,xlon(100),xlon(110),xlat(10),xlat(20)
        Deallocate (xtmp,lon,lat,xt,XX,xave) !,xtmp1,xtmp2,ymean)
      End Subroutine Read_var
      !-------------------------------------------------------------------------------------------
      Subroutine Create_file_for_Statistical_Models (ilt,Out_Dir) !(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)

         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, ilt
         Character (Len=7) :: EOFName(1:NMon,1:Nmode)
         Character (Len=255) :: Fname, Fname1, Fname2
         Character (Len=10)  :: vnametc
         Character (Len=*) :: Out_Dir

         vnametc=trim(adjustl(vartc))

         Do ieigen=1, Nmode
            Do imon=1, NMon
               write(EOFName(imon,ieigen),"(A3,I1,A1,I2.2)") "EOF",ieigen,"_",imon
            Enddo
         Enddo

         Do imon=1, NMon
            !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 ", imon, " and Predictors"
            write (3,"(3I6)") Yr1+1, Nt, 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,Nt    ! EOFall(LT,Nmon,Nmode,Nt)
               write (3,"(20F10.4)") TC(iyr+1,imon), EOFall(ilt,imon,iyr,1:Nmode)
            Enddo
            Close (3)
         Enddo
         !Close (33); Close (44); Close (55)
      End Subroutine Create_file_for_Statistical_Models
      !---------------------------------------------------------------------------------
      Subroutine Calc_EOFs_for_Fcs_DL !(Fname)
        !Character (Len=*) :: Fname
        Real, Allocatable :: lon(:),lat(:),lev(:),xtime(:),xtmp(:)
        Real, Allocatable :: Xpre(:), Fpre(:), A(:,:), x11(:), Fx(:,:), Tmp
        Real, Allocatable :: x2d(:,:), y2d(:,:), lamda(:,:), tb(:),sx(:)
        Integer :: nlon,nlat,nz,ntime
        Integer :: ilt,imon,it,imod

        it1=20     ! Lay tu 2003, nhay qua 20 nam dau, tu 1983-2002
        Nt1=Nt-it1 ! Lay tu 2003, nhay qua 20 nam dau, tu 1983-2002
        Allocate (xtmp(Nx),x11(Nmode))
        Allocate (Xpre(Nx), Fpre(Nmode), A(Nx,Nmode))
        Allocate (EOFall(LT,Nmon,Nt1,Nmode))

        Allocate (x2d(Nt1,Nx),y2d(Nt1,Nx),lamda(Nmode,Nmode),tb(Nx),sx(Nx),Fx(Nt1,Nmode))

        Do ilt=1,LT
           Do imon=1,Nmon
              tmp=Sqrt(Nsple(ilt,imon,1))
              !A(1:Nx,1:Nmode)=PCall(ilt,imon,1:Nx,1:Nmode)    ! Cac PCs
              Do j=1,Nx; A(j,1:Nmode)=PCall(ilt,imon,j,1:Nmode); Enddo  ! A(Nx,K)
              Do j=1,Nx; tb(j)=xtball(ilt,imon,j); sx(j)=sxall(ilt,imon,j); Enddo

              Do i=1,Nt1; x2d(i,1:Nx) = Xfc(ilt,imon,i+it1,1:Nx); Enddo

              Do it=1,Nt1
                 Do j=1,Nx
                    y2d(it,j) = (x2d(it,j)-tb(j))/sx(j)/tmp    ! Y(Nt,Nx)=(X(Nt,Nx)-Xtb(Nx))/Sx(Nx)/sqrt(Nt)
                 Enddo
              Enddo
              lamda = 0
              Do j=1,Nmode
                 lamda(j,j) = 1./Eigall(ilt,imon,j)            ! 1/sqrt(eig)
              Enddo

              Fx = Matmul(Matmul(y2d,A),lamda)                 ! F(Nt,K) = Y(Nt,Nx) x A(Nx,K) x lamda(K,K)
              Do i=1,Nt1
                 EOFall(ilt,imon,i,1:Nmode)=Fx(i,1:Nmode)
              Enddo
           Enddo
        Enddo
        !print*,EOFall
        print*,"Calc_EOFs_for_Fcs_DL...Done"
      End Subroutine Calc_EOFs_for_Fcs_DL

      !-------------------------------------------------------------------------------------------
      Subroutine Create_file_for_Fcs_DLap (ilt,Out_Dir) !(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)

         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, ilt
         Character (Len=7) :: EOFName(1:NMon,1:Nmode)
         Character (Len=255) :: Fname, Fname1, Fname2
         Character (Len=10)  :: vnametc
         Character (Len=*) :: Out_Dir

         vnametc=trim(adjustl(vartc))

         print*,"In Create Files: Nmode,NMon,Nt= ",Nmode,NMon,Nt
         Do ieigen=1, Nmode
            Do imon=1, NMon
               write(EOFName(imon,ieigen),"(A3,I1,A1,I2.2)") "EOF",ieigen,"_",imon
            Enddo
         Enddo

         Do imon=1, NMon
            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")

            it1=20;
            Nt1=Nt-it1 ! Tinh tu 2003-2009
            write (3,"(A,I2.2,A)") "Total TC count of month ", imon, " and Predictors"
            write (3,"(3I6)") Yr1+it1+1, Nt1, Nmode+1
            write (3,"(20I3)") 3, (0,ivar=1,Nmode)   !(IDX(i), i=1,Mvar)
            write (3,"(20A10)") trim(vnametc),(trim(EOFName(imon,i)),i=1,Nmode)
            Do iyr=1,Nt1    ! EOFall(LT,Nmon,Nmode,Nt)
               write (3,"(20F10.4)") TC(iyr+it1+1,imon), EOFall(ilt,imon,iyr,1:Nmode)
            Enddo
            Close (3)
         Enddo
         !Close (33); Close (44); Close (55)
      End Subroutine Create_file_for_Fcs_DLap
      !-------------------------------------------------------------------------------------------

      Include "/home/tanpv/Mylibs/Sub_PCA_2019.f90"
      Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
      Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
      !Include "Write_EOF_PCA_for_Plot.f90"
      !Include "Write_EOF_PCA_for_Fcs.f90"
      !Include "Create_file_for_Statistical_Models.f90"
      !-------------------------------------------------------------------------------------------
      END

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