PROGRAM xevlmem
      Use NetCDF
      Real, Parameter :: Miss=-99.0
      Integer, Parameter :: YrBL1=1981, YrBL2=2010, Nmon=12
      ! While exact values depend on the application, 
      ! one might take Mp = 10 or 20 or 50 for Ntime = 1000 or 10000
      !Integer, Parameter :: Mp = 20

      Character (Len=150) :: Fname
      REAL,ALLOCATABLE :: Spec(:), Freq(:), Peri(:)
      Real, Allocatable :: TC_NWP(:), Ann_Mean(:), Std(:)

      Real, Allocatable :: cof(:)

      !INTEGER :: N, NFDT, M
      !PARAMETER(N=1000,M=100,NFDT=16)

      INTEGER :: i
      REAL :: fdt, pm

      Integer :: Ntime, Nf, Mp, Time_Step
     !----------------------------------------------------------------------------------------------------

      Call Read_TC_NC
      Time_Step=3; Mp=50
      Nf=Ntime/Time_Step
      
      Allocate (cof(Mp))
      Allocate (Spec(Nf), Freq(Nf), Peri(Nf))

      !call memcof(data,N,M,pm,cof)
      call memcof(TC_NWP,Ntime,Mp,pm,cof)
      write(*,'(1x,t6,a,t20,a)') 'f*delta','power'
      do i=1, Nf !NFDT
        fdt=0.5*i/Nf !NFDT
        Freq(i)=fdt; Peri(i)=1.0/fdt*Time_Step/Nmon
        Spec(i)=evlmem(fdt,cof,Mp,pm)
        write(*,'(1x,i4,3f12.6)') i, log(Peri(i)),Freq(i),Spec(i) !evlmem(fdt,cof,Mp,pm)
      enddo
      Call Write_Spec_TC_to_NC

     CONTAINS 
     !----------------------------------------------------------------------------------------------------
        include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
        include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"

     !------------------------------------------------------------------------------
      Subroutine Write_Spec_TC_to_NC
        Integer :: iw, ifre
        Real, Allocatable :: PSD(:,:), xlat(:), Xfre(:), Per(:), xlon(:)

        Allocate (PSD(Nf,1), xlat(1), Per(Nf)) !Xfre(1:MP1),Peri(1:MP1), xlon(MP1))
        Do ifre=1, Nf
           PSD(ifre,1)=Spec(Nf-ifre+1)
           Per(ifre)  =Peri(Nf-ifre+1)
        Enddo
        !Do ifre=1, MP1
        !   Peri(MP1-ifre+1)=Log(Period(ifre))
        !   Xfre(ifre) = Freq(ifre)
        !   xlon(ifre) = ifre
        !Enddo
        Fname="Spectral_TC.nc"
        xlat(1)=1.0
        Call wr_ncdf_2d (Fname, Per, xlat, Nf, 1, PSD, "spec", &
                         "spectral density", "-", "peri", "lat", Miss)
!        Call wr_ncdf_2d (Fname, Freq, xlat, Nf, 1, PSD, "spec", &
!                         "spectral density", "-", "fre", "lat", Miss)
        !Call wr_ncdf_1d_apd (trim(Fname), Peri, "peri", Nf, &
        !     "Period", "mon", "fre", Miss)
        !Call wr_ncdf_1d_apd (trim(Fname), Peri, "peri", MP1, &
        !     "period", "-", "stt", Miss)

      End Subroutine Write_Spec_TC_to_NC
     !----------------------------------------------------------------------------------------------------
      Subroutine Read_TC_NC
        Integer :: Nyear, Nmon, iyr, imon, it
        Real, Allocatable :: TC(:,:), Year(:), Mon(:)
        Character (Len=150) :: Fname
        Integer, Parameter :: Mave=6   !  Step of Moving summation
        Real :: tmp

        !Fname="/home/tanpv/OBS/KC09/TC_mon.nc"
        Fname="/home/tanpv/OBS/KC09/VES_Monthly_TC_NTC_Count_jp.nc"
        !Call rd_ncdf_2d (Fname,Year,Mon,Nyear,Nmon, TC, "jma_ves_tc") !"JMA_NWP_TC_All")
        Call rd_ncdf_2d (Fname,Year,Mon,Nyear,Nmon, TC, "tc")
        Ntime=Nyear*Nmon
        Allocate (TC_NWP(Ntime), Ann_Mean(Nmon), Std(Nmon))
        ! Climatological Annual Cycle:
        Ann_Mean(:)=0; Std(:)=0
        Do imon=1, Nmon
           Do iyr=YrBL1, YrBL2
              it=iyr-Year(1)+1; !print*, it,TC(it,imon)
              Ann_Mean(imon)=Ann_Mean(imon)+TC(it,imon)
              !Std(imon)=TC(it,imon)*TC(it,imon)
           Enddo
           Ann_Mean(imon)=Ann_Mean(imon)/(YrBL2-YrBL1+1)
           !Std(imon)=SQRT(Std(imon)/(YrBL2-YrBL1+1)-Ann_Mean(imon)*Ann_Mean(imon))
        Enddo
        ! Anomalous timeseries:
        it=0
        Do iyr=1, Nyear
           Do imon=1, Nmon
              it=it+1; TC_NWP(it)=(TC(iyr,imon)-Ann_Mean(imon)) !/Std(imon)
           Enddo
        Enddo
        print*,"TC data: ",Nyear,Nmon, Ntime
      End Subroutine Read_TC_NC
     !------------------------------------------------------------------------------

      SUBROUTINE memcof(data,n,m,xms,d)
      INTEGER :: m,n !,MMAX,NMAX
      REAL :: xms,d(m),data(n)
      !PARAMETER (MMAX=60,NMAX=2000)
      INTEGER :: i,j,k
      REAL :: denom,p,pneum !,wk1(NMAX),wk2(NMAX),wkm(MMAX)
      Real :: wk1(n), wk2(n), wkm(m)
      !if (m.gt.MMAX.or.n.gt.NMAX) pause 'workspace too small in memcof'
      p=0.
      do j=1,n
         p=p+data(j)**2
      enddo
      xms=p/n
      wk1(1)=data(1)
      wk2(n-1)=data(n)
      do j=2,n-1
         wk1(j)=data(j)
         wk2(j-1)=data(j)
      enddo
      do k=1,m
         pneum=0.
         denom=0.
         do j=1,n-k
            pneum=pneum+wk1(j)*wk2(j)
            denom=denom+wk1(j)**2+wk2(j)**2
         enddo
         d(k)=2.*pneum/denom
         xms=xms*(1.-d(k)**2)
         do i=1,k-1
            d(i)=wkm(i)-d(k)*wkm(k-i)
         enddo
         if(k.eq.m)return
         do i=1,k
            wkm(i)=d(i)
         enddo
         do j=1,n-k-1
            wk1(j)=wk1(j)-wkm(k)*wk2(j)
            wk2(j)=wk2(j+1)-wkm(k)*wk1(j+1)
         enddo
      enddo
      pause 'never get here in memcof'
      END SUBROUTINE memcof
      !------------------------------------------------------------------------
      Real FUNCTION evlmem(fdt,d,m,xms)
      INTEGER :: m
      REAL :: fdt,xms,d(m)
      INTEGER :: i
      REAL :: sumi,sumr
      DOUBLE PRECISION :: theta,wi,wpi,wpr,wr,wtemp

      theta=6.28318530717959d0*fdt
      wpr=cos(theta)
      wpi=sin(theta)
      wr=1.d0
      wi=0.d0
      sumr=1.
      sumi=0.
      do i=1,m
        wtemp=wr
        wr=wr*wpr-wi*wpi
        wi=wi*wpr+wtemp*wpi
        sumr=sumr-d(i)*sngl(wr)
        sumi=sumi-d(i)*sngl(wi)
      enddo
      evlmem=xms/(sumr**2+sumi**2)
      return
      END FUNCTION evlmem
      !------------------------------------------------------------------------
      END Program