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, vartc
      Character (Len=1)   :: Lst
      !-------------------------------------------------------------------------------------------
      Real, Allocatable :: Corre(:,:,:,:), Pval(:,:,:,:) ! (Nlon,Nlat,Nmon,LT)
      Real, Allocatable :: X1D(:), X2D(:,:), X3D(:,:,:)  ! (Nlon,Nlat,Nt)
      !-------------------------------------------------------------------------------------------
      Integer, Parameter :: Nvar=7
      Integer :: lev, ivar, Mov_Ave
      ! AO fields:
      Character (Len=5) :: Vars(Nvar)= &
               (/"sst","prmsl","vor","vws","vvel","rh","h85"/)
      !                                     850,   700, 850
      Character (Len=20), Allocatable :: Varn_out(:,:,:)
      Integer :: Nvung(Nvar)
      Integer, Allocatable :: Targ_mon(:), Mvar(:,:), Mvung(:,:,:), LeadTime(:)
      Real, Allocatable :: Vung1(:,:,:,:), Vung2(:,:,:,:), Xfc(:,:,:,:)  ! So lieu cho mo hinh du bao
      Integer :: ii, jj, kk, Mvarx
      !-------------------------------------------------------------------------------------------
      
      Mvarx=Nvar*2 + 1    ! 7 vars x 2 vung + TC

      Open (1,file="Extract_data.txt", status="old")
      read (1,*)
      Read (1,*) NLT; Allocate (Leadtime(NLT))
      Do jlt=1, NLT
         read (1,*) ilt; Leadtime(jlt)=ilt
         read (1,*) Ntag_mon     ! So thang lam du bao (Max=Nmon)
         Allocate (Mvar(NLT,Ntag_mon), Xfc(NLT,Ntag_mon,Nyrs,Mvarx))
         Allocate (Mvung(NLT,Ntag_mon,Nvar))
         Allocate (Varn_out(NLT,Ntag_mon,Mvarx)); Varn_out=""
         print*, NLT, ilt, Ntag_mon
         Allocate (Vung1(NLT,Ntag_mon,Nvar,4),Vung2(NLT,Ntag_mon,Nvar,4))
         Allocate (Targ_mon(Ntag_mon))
         Do imon=1, Ntag_mon
            Mvar(jlt,imon)=0; kk=0
            read (1,*) Targ_mon(imon)
            print*,Targ_mon(imon)
            Do ivar=1,Nvar
               read (1,*) ii, iv; Mvung(jlt,imon,ivar)=iv
               if (iv==1) then
                  kk=kk+1; Varn_out(jlt,imon,kk)=trim(Vars(ivar))//"_v1"
               else if (iv==2) then
                  kk=kk+1; Varn_out(jlt,imon,kk)=trim(Vars(ivar))//"_v1"
                  kk=kk+1; Varn_out(jlt,imon,kk)=trim(Vars(ivar))//"_v2"
               endif
               Mvar(jlt,imon)=Mvar(jlt,imon)+iv    ! So bien doc lap = Nvar * So vung
               !print*, ii, iv; pause
               backspace(1)
               read (1,*) ii, jj, (Vung1(jlt,imon,ivar,lonlat),lonlat=1,4), &
                                  (Vung2(jlt,imon,ivar,lonlat),lonlat=1,4)
               !print"(3I3,10F6.0)",jlt, ivar, iv, Vung1(jlt,imon,ivar,:),Vung2(jlt,imon,ivar,:)
            Enddo
            !print"(14A10)",Varn_out(jlt,imon,:)
         Enddo
      Enddo
      Close (1)

      Fname="/home/tanpv/Data/TCs/VES_Monthly_TC_NTC_Count.nc"
      vartc="tc3m"
      Mov_Ave=3     !  =0 neu khong tinh trung binh truot (cho monthly)
                    ! /=0 Tinh trung binh mua (3 hoac 6 thang)
      ! Reading TC data:
      Call Read_TC (Fname,trim(vartc))
      
      ! Reading AO fields:
      Dir="/home/tanpv/Data/CFS_Rfc/"//trim(Subdir) !ensmean_all/"
      !---------------------------------------------------------------------------------
      Do ivar=1,Nvar
         vname=trim(Vars(ivar))
         print*,"Calculating for ",vname
         if (ivar <= 4) then      ! Reading 3D-arrays (lon,lat,time)
            Call Read_var (vname,Mov_Ave)
         else if (ivar<=7) then   ! Reading 4D-arrays (lon,lat,lev,time)
            if (trim(Vars(ivar))=="h85") then
               lev=850
               Call Read_var_4D ("hgt",lev,Mov_Ave)
            else if (trim(Vars(ivar))=="vvel") then
               lev=850
               Call Read_var_4D (vname,lev,Mov_Ave)
            else if (trim(Vars(ivar))=="rh") then
               lev=700
               Call Read_var_4D (vname,lev,Mov_Ave)
            endif
            !Call Read_var_4D (vname,lev,Mov_Ave)
         !else if (ivar==7) then        ! Reading 4D-arrays, multi levs (700,600,500)
         !   Call Read_var_4D_mlv ("rh",Mov_Ave)
         !else if (ivar==9) then   ! h850
         !   lev=850
         !   Call Read_var_4D ("hgt",lev,Mov_Ave)
         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
      Call Extract_Var
      Call Write_to_TXT

      !---------------------------------------------------------------------------------
      CONTAINS
      !-------------------------------------------------------------------------------------------
      Subroutine Write_to_TXT
        Integer :: it,ilt,imon,ilon,ilat,ix,iy,ivar,iv,jlt,mva
        Character (Len=100) :: Fname, Dir
        ! Varn_out(jlt,imon,kk)
        ! Xfc(jlt,it,1:Nyrs,iv)

        Dir="/home/tanpv/Progs/TCs/FCS/Dat_out/"
        Do jlt=1,NLT
           Do it=1,Ntag_mon
              imon=Targ_mon(it)
              write (Fname,"(A,I1,A,I2.2,A)") trim(Dir)//trim(vartc)// &
                           "_LT_",Leadtime(jlt),"_Mon_",imon,".txt"
              Open (3,file=Fname,status="Unknown")
              Mv=Mvar(jlt,it)
              write (3,"(2A,I4,A,I4)") "So lieu de du bao "//trim(vartc), &
                          " thoi ky ",Yr1,"-",Yr2
              write (3,"(A,I2.2,A,I1,A)") "(Thang du bao: ",imon, &
                                ". Han du bao: ",Leadtime(jlt)," thang)"
              write (3,"(3A)") "(Cot 1: TC, Cac cot con lai la cac bien", &
                               "SST:C; Pmsl:mb; vor:*10^6; vws:m/s", &
                               " vvel:*10^2; rh:%; h85:dam; v1/v2:Vung)"
              write (3,"(3I6)") Yr1, Nyrs, Mv+1
              write (3,"(20I10)") 3, (ix-ix, ix=1,Mv)
              write (3,"(20A10)") trim(vartc), (trim(Varn_out(jlt,it,iv)), iv=1,Mv)
              Do iyr=1,Nyrs
                 write (3,"(20F10.4)") (Xfc(jlt,it,iyr,iv),iv=1,Mv+1)  ! Ca TC
              Enddo
              Close (3)
           Enddo
        Enddo
        ! Write out for PCA:
        Do jlt=1,NLT
           write (Fname,"(A,I1,A)") trim(Dir)//"TC_PCA_Leadtine_",Leadtime(jlt),".txt"
           Open (3,file=Fname, status="Unknown")
        !   write (3,"(I3,A)") NLT,"    ! So Leadtime"
           write (3,"(I3,A)") Leadtime(jlt), "       ! Leadtime" 
           write (3,"(I6,A)") Ntag_mon,"      ! So thang du bao"      
           Do it=1,Ntag_mon
              imon=Targ_mon(it)
              Mv=Mvar(jlt,it)
              !write (3,"(I6,A)") Ntag_mon,"      ! So thang du bao"
              write (3,"(I6,A)") imon,    "      ! Thang duoc du bao"
              write (3,"(2I6,A)") Nyrs, Mv+1,"   ! Dg lg mau, so bien"
              write (3,"(20A10)") trim(vartc), (trim(Varn_out(jlt,it,iv)), iv=1,Mv)
              Do iyr=1,Nyrs
                 write (3,"(20F10.4)") (Xfc(jlt,it,iyr,iv),iv=1,Mv+1)  ! Ca TC
              Enddo
           Enddo
           Close (3)
        Enddo
      End Subroutine Write_to_TXT
      !-------------------------------------------------------------------------------------------
      Subroutine Extract_Var
        Integer :: it,ilt,imon,ilon,ilat,ix,iy,ivar,iv,jlt
        Real :: xlon1,xlon2,xlat1,xlat2

        !Xfc(NLT,Ntag_mon,Nyrs,Mvarx), Targ_mon(Ntag_mon)
        !Xall(:,:,:,:,:,:) ! (Nlon,Nlat,Nyrs,Nmon,LT,Nvar)
        Allocate (X1D(Nyrs),X2D(Nlon,Nlat),X3D(Nlon,Nlat,Nyrs))

        Do jlt=1,NLT
           ilt=Leadtime(jlt)
           Do it=1,Ntag_mon
              imon=Targ_mon(it); iv=1    ! iv=1 for TC data
              Do ivar=1,Nvar
                 X3D(1:Nlon,1:Nlat,1:Nyrs)=Xall(1:Nlon,1:Nlat,1:Nyrs,imon,ilt,ivar)
                 !print*,X3D(130,110,3:5);pause
                 if (Mvung(jlt,it,ivar)==1) then         ! Chi lay 1 vung
                    xlon1=Vung1(jlt,it,ivar,1); xlon2=Vung1(jlt,it,ivar,2)
                    xlat1=Vung1(jlt,it,ivar,3); xlat2=Vung1(jlt,it,ivar,4)
                    Call Subdomain (xlon1,xlon2,xlat1,xlat2)  ! Chuoi Trung binh vung 1: X1D
                    iv=iv+1; Xfc(jlt,it,1:Nyrs,iv)=X1D(1:Nyrs)
                 else if (Mvung(jlt,it,ivar)==2) then    ! Lay 2 vung
                    xlon1=Vung1(jlt,it,ivar,1); xlon2=Vung1(jlt,it,ivar,2)
                    xlat1=Vung1(jlt,it,ivar,3); xlat2=Vung1(jlt,it,ivar,4)
                    Call Subdomain (xlon1,xlon2,xlat1,xlat2)  ! Chuoi Trung binh vung 1: X1D
                    iv=iv+1; Xfc(jlt,it,1:Nyrs,iv)=X1D(1:Nyrs)

                    xlon1=Vung2(jlt,it,ivar,1); xlon2=Vung2(jlt,it,ivar,2)
                    xlat1=Vung2(jlt,it,ivar,3); xlat2=Vung2(jlt,it,ivar,4)
                    Call Subdomain (xlon1,xlon2,xlat1,xlat2)  ! Chuoi Trung binh vung 1: X1D
                    iv=iv+1; Xfc(jlt,it,1:Nyrs,iv)=X1D(1:Nyrs)
                 endif
                 !print*,"imon=",imon," ivar=",ivar; pause
              Enddo
           Enddo
        Enddo
      End Subroutine Extract_Var
      !-------------------------------------------------------------------------------------------
      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,Mov_Ave)
        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
        Integer :: Mov_Ave

        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
           if (Mov_Ave>0) then  ! Trung binh truot
              Do ilat=1,Nlat
                 Do ilon=1,Nlon
                    xave(1:Ntime)=xtmp(ilon,ilat,1:Ntime)
                    Call Moving_Ave (xave,Ntime,Mov_Ave)
                    xtmp(ilon,ilat,1:Ntime)=xave(1:Ntime)
                 Enddo
              Enddo
           endif
           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,Mov_Ave)
        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
        Integer :: Mov_Ave

        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
           if (Mov_Ave>0) then  ! Trung binh truot
              Do ilat=1,Nlat
                 Do ilon=1,Nlon
                    xave(1:Ntime)=xtmp(ilon,ilat,1:Ntime)
                    Call Moving_Ave (xave,Ntime,Mov_Ave)
                    xtmp(ilon,ilat,1:Ntime)=xave(1:Ntime)
                 Enddo
              Enddo
           endif
           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
      !-------------------------------------------------------------------------------------------
      Subroutine Moving_Ave (X,N,m)
         ! Moving average with m steps and filling up all first X's elements
         Integer :: N,m,i,j,mj
         Real :: X(N), Y(N), tmp

         Y(:)=Miss
         Do i=1,N-m+1
            tmp=0; mj=0
            Do j=i,i+m-1
               if (X(j)/=Miss) then
                  mj=mj+1; tmp=tmp+X(j)
               endif
            Enddo
            if (mj==m) Y(i)=tmp/mj
         Enddo
         X(:)=Y(:)
      End Subroutine Moving_Ave
      !-------------------------------------------------------------------------------------------
      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,Mov_Ave)
        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(:)
        Integer :: Mov_Ave
        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(:)
              !Nlon=Nx; Nlat=Ny;
              !Allocate (xlon(Nlon),xlat(Nlat),xtmp(Nlon,Nlat,Ntime))
              !Allocate (sst(Nlon,Nlat,Nyrs,Nmon,LT),xave(Ntime))

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

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

        !Nlon=Int((Lon2-Lon1)/Dlon)+1; Nlat=Int((lat2-Lat1)/Dlat)+1
        !ix=Int(Lon1/Dlon)+1; iy=Int((Lat1+90)/Dlat)+1
        !Allocate (xlon(Nlon),xlat(Nlat),xtmp(Nlon,Nlat,Ntime))
        !Allocate (sst(Nlon,Nlat,Nyrs,Nmon,LT),xave(Ntime))
        !Allocate (X2D(Nlon,Nlat))

        !xlon(1:Nlon) = lon(ix:ix+Nlon); xlat(1:Nlat) = lat(iy:iy+Nlat)
        ! Calculate the moving average first:

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

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

        !Allocate (xtmp(1:Nx,1:Ny,1:Ntime),xave(Ntime))
        !Allocate (sst(Nx,Ny,Nyrs,Nmon,LT))
        !Allocate (X2D(1:Nx,1:Ny), X3D(1:Nx,1:Ny,1:Nt), X1D(Nt))
        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
           if (Mov_Ave>0) then  ! Trung binh truot
              Do ilat=1,Nlat
                 Do ilon=1,Nlon
                    xave(1:Ntime)=xtmp(ilon,ilat,1:Ntime)
                    Call Moving_Ave (xave,Ntime,Mov_Ave)
                    xtmp(ilon,ilat,1:Ntime)=xave(1:Ntime)
                 Enddo
              Enddo
           endif
           !
           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

        ! Test with SST:
        if (trim(vname)=="sst") then
           Call Gradient
        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 Subdomain (xlon1,xlon2,xlat1,xlat2)
         Real :: xmean,xlon1,xlon2,xlat1,xlat2, tmp
         Integer :: iyr, ilon,ilat,ic, Nx,Ny

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

         X1D(:)=Miss
         Do iyr=1, Nyrs
            X2D(:,:)=X3D(:,:,iyr)
            tmp=0.0; ic=0
            Do ilon=ix, ix+Nx
               Do ilat=iy, iy+Ny
                  if (X2D(ilon,ilat) /= Miss) then
                     ic=ic+1; tmp=tmp+X2D(ilon,ilat)
                  endif
               Enddo
            Enddo
            if (ic>0) then
               xmean=tmp/real(ic)
            else
               xmean=Miss
            endif
            X1D(iyr)=xmean  ! Chuoi thoi gian trung binh mien
            !print*,xlon1,xlon2,xlat1,xlat2
         Enddo
      End Subroutine Subdomain
      !-------------------------------------------------------------------------------------------
      Subroutine Gradient
         Integer :: ilon,ilat,iyr,imon,ilt
         Real, Allocatable :: xtmp1(:,:,:,:,:),xtmp2(:,:,:,:,:)
         Real :: tmp

         Allocate (xtmp1(Nlon,Nlat,Nyrs,Nmon,LT))
         Allocate (xtmp2(Nlon,Nlat,Nyrs,Nmon,LT))
         Allocate (X2D(Nlon,Nlat))
         xtmp1=sst; xtmp2=Miss

         Do ilt=1,LT
            Do iyr=1,Nyrs
               Do imon=1,Nmon
                  Do ilon=1,Nlon; Do ilat=1,Nlat
                     X2D(ilon,ilat)=xtmp1(ilon,ilat,iyr,imon,ilt)
                  Enddo; Enddo
                  tmp=Sub_mean (imon)
                  !tmp=xtmp1(69,34,iyr,imon,ilt)
                  if (tmp /= Miss) then
                     Do ilon=1,Nlon; Do ilat=1,Nlat
                        if (xtmp1(ilon,ilat,iyr,imon,ilt)/=Miss) &
                            xtmp2(ilon,ilat,iyr,imon,ilt) = &
                            xtmp1(ilon,ilat,iyr,imon,ilt)-tmp
                     Enddo; Enddo
                  endif
               Enddo
            Enddo
         Enddo
         sst=xtmp2
         Deallocate (xtmp1,xtmp2,X2D)
      End Subroutine Gradient
      !-------------------------------------------------------------------------------------------
      Real Function Sub_mean (imon)
        Character (Len=3) :: SubDomain
        Integer :: ilon,ilat, ilo1,ilo2,ila1,ila2
        Integer :: imon,iyr,ic, ix, iy
        Real :: tmp, xmean, xlon1,xlon2,xlat1,xlat2

        !      Lo1-Lo2  La1-La2  i1-i2   j1-j2
        ! PLP: 130-140,   8- 18 (71-81,  29-39)
        ! VES: 110-118,  10- 18 (51-59,  31-39)
        ! CHN: 130-140,  20- 30 (71-81,  41-51)
        ! IND:  80- 90,  -5-  5 (21-31,  16- 26)
        !ilo1=70; ilo2=75; ila1=30; ila2=35

        Select Case (imon)
           Case (1,6,7)
              SubDomain="PLP"
              !ilo1=71; ilo2=81; ila1=29; ila2=39
              xlon1=128; xlon2=129; xlat1=9; xlat2=10
              Nx=Int((xlon2-xlon1)/Dlon)+1; Ny=Int((xlat2-xlat1)/Dlat)+1
              ix=Int(xlon1/Dlon)+1; iy=Int((xlat1+90)/Dlat)+1
              ilo1=ix; ilo2=ix+1; ila1=iy; ila2=iy+1
           Case (2,3,4,5,11,12)
              SubDomain="PLP" !"VES"
              !ilo1=51; ilo2=59; ila1=31; ila2=39
              xlon1=115; xlon2=116; xlat1=14; xlat2=15
              Nx=Int((xlon2-xlon1)/Dlon)+1; Ny=Int((xlat2-xlat1)/Dlat)+1
              ix=Int(xlon1/Dlon)+1; iy=Int((xlat1+90)/Dlat)+1

              ilo1=ix; ilo2=ix+1; ila1=iy; ila2=iy+1
           Case (8)
              SubDomain="PLP" !"CHN"
              !ilo1=71; ilo2=81; ila1=41; ila2=51
              !ilo1=68; ilo2=70; ila1=41; ila2=43
              xlon1=85; xlon2=86; xlat1=-5; xlat2=-4
              Nx=Int((xlon2-xlon1)/Dlon)+1; Ny=Int((xlat2-xlat1)/Dlat)+1
              ix=Int(xlon1/Dlon)+1; iy=Int((xlat1+90)/Dlat)+1

              ilo1=ix; ilo2=ix+1; ila1=iy; ila2=iy+1  ! IND
           Case (9,10)
              SubDomain="PLP" !"IND"
              !ilo1=21; ilo2=31; ila1=16; ila2=26
              xlon1=85; xlon2=86; xlat1=-5; xlat2=-4
              Nx=Int((xlon2-xlon1)/Dlon)+1; Ny=Int((xlat2-xlat1)/Dlat)+1
              ix=Int(xlon1/Dlon)+1; iy=Int((xlat1+90)/Dlat)+1

              ilo1=ix; ilo2=ix+1; ila1=iy; ila2=iy+1
           Case DEFAULT
              SubDomain="XXX"
              Sub_mean=0; Return
        End Select

        tmp=0.0; ic=0
        Do ilat=ila1,ila2
           Do ilon=ilo1,ilo2
              if (X2D(ilon,ilat) /= Miss) then
                 ic=ic+1; tmp=tmp+X2D(ilon,ilat)
              endif
           Enddo
        Enddo
        if (ic>0) then
           xmean=tmp/real(ic)
        else
           xmean=Miss
        endif
        !print*,ilo1,ilo2,ila1,ila2,xmean; pause
        Sub_mean=xmean
      End Function Sub_mean

      !-------------------------------------------------------------------------------------------
      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 "/home/tanpv/Mylibs/Calc_Correlation_Coefficient.f90"
      !-------------------------------------------------------------------------------------------
      END
      Include "/home/tanpv/Mylibs/Student_Prob.f90"
      !-------------------------------------------------------------------------------------------