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

        print*,"Begining Calc_EOFs_for_Fcs_DL"
        Call rd_ncdf_4d (Fname,lon,lat,lev,xtime,nlon,nlat,nz,ntime,PCall,"pc")
        Nmode=ntime
        Deallocate (lon,lat,lev,xtime)
        Call rd_ncdf_3d (Fname,lon,lat,lev,nlon,nlat,nz,Eigall,"eig")
        Deallocate (lon,lat,lev)
        Call rd_ncdf_3d (Fname,lon,lat,lev,nlon,nlat,nz,xtball,"mean")
        Deallocate (lon,lat,lev)
        Call rd_ncdf_3d (Fname,lon,lat,lev,nlon,nlat,nz,sxall,"std")
        Deallocate (lon,lat,lev)
        Call rd_ncdf_3d (Fname,lon,lat,lev,nlon,nlat,nz,Nsple,"nt")
        Deallocate (lon,lat,lev)
        print*,"Reading data done !!!"
        print*, "Nx,Nt,Nmode,LT,Nmon=",Nz,Nt,Nmode,LT,Nmon

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