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 :: Nmon=12, Yr1=2011, Yr2=2017, 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=".2011-2017."
      Character (Len=15), Parameter :: Subdir="cfs1117/"

!      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(:,:), TC1(:,:)        ! (Nyrs,Nmon)
      Real, Allocatable :: xlon(:),xlat(:),xyrs(:),xmon(:)  ! Nlon,Nlat,Nyrs,Nmon
      Character (Len=100) :: Fname, Dir
      Character (Len=10)  :: vname, vartc
      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=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
      !Namelist /Model_Info/ Fname, vname, pca_out, eof_out, &
      !                      Yr1, Yr2, Nmon, &
      !Integer :: switch, standard, PC_Stard, Nmode
      
      !-------------------------------------------------------------------------------------------
      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="tc6m"
      ! Reading TC data:
      Call Read_TC (Fname,trim(vartc))
      
      ! Reading AO fields:
      Dir="/home/tanpv/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
      Enddo
      print*,"Reading all data ... Done !"
      Call Extract_Var_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      =  5         ! So mode duoc giu lai (Nmode <= Min{nt,nx})

      Fname="Dat_out/pca_tc3m_all.nc" 
      Call Calc_EOFs_for_Fcs_DL (Fname)

      ilt=3
      Out_Dir="/home/tanpv/Progs/TCs/FCS/Out_Fcs/"
      Call Create_file_for_Fcs_DLap (ilt,Out_Dir)
      !---------------------------------------------------------------------------------
      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+1 -> 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
      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

        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     ! 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

        ! level = 850, 700, 500, 300, 200 ;

        if (lev==850) then
           ilev=1
        elseif (lev==700) then
           ilev=2
        elseif (lev==500) then
           ilev=3
        elseif (lev==300) then
           ilev=4
        elseif (lev==200) then
           ilev=5
        endif
        Level=ilev
      End Function Level
      !-------------------------------------------------------------------------------------------
      Integer Function Level_82_09 (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_82_09=ilev
      End Function Level_82_09
      !-------------------------------------------------------------------------------------------
      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
      !-------------------------------------------------------------------------------------------
      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 "Calc_EOFs_for_Fcs_DL_New.f90"
      Include "Create_file_for_Fcs_DLap.f90"
      !-------------------------------------------------------------------------------------------
      END

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