! CHUONG TRINH NAY TINH PHO ENTROPI CUC DAI 
! CHO SO LIEU CACH NHAU MOT NAM (IYEAR=1) HOAC SL THANG (IYEAR=0)
! MODIFY ON 6/5/2000
      PROGRAM MEMSP
      Use NetCDF
      Real, Parameter :: Miss=-99.0
      REAL PM
      REAL DELTA_S
      INTEGER M, MP, MP1
!
      Character (Len=150) :: Fname
      REAL,ALLOCATABLE :: DAT1(:),COF(:),S(:),FREQ(:)
      REAL,ALLOCATABLE :: WK1(:),WK2(:),WKM(:)
      REAL,ALLOCATABLE :: DAT(:,:)       ! MANG CHUA SO LIEU BAN DAU
      Integer :: Ntime, Nwave            ! To be provided by Read_Data Sub.
      Real, Allocatable :: imf(:,:), Coff(:,:), Sf(:,:)  !, Freqf(:,:)
      Real, Allocatable :: SMean(:), Dfre(:), Period(:)
      Integer :: I, iw
! ----------------------------------------------------------------------------------
      !DELTA_S = 0.001
      !MP  = 0.5/DELTA_S   ! fNyq: Nyquist frequency
      !MP1 = MP+1

      M   = 30          ! Need to be tested
! ----------------------------------------------------------------------------------
      Call Read_Data  ! Return Ntime, Nwave, imf(Ntime,Nwave)
      DELTA_S = 1./Real(Ntime)
      MP  = 0.5/DELTA_S
      MP1 = MP+1
! ----------------------------------------------------------------------------------
      ALLOCATE (DAT1(Ntime),COF(M),S(MP1),FREQ(MP1))
      !ALLOCATE (WK1(Ntime),WK2(Ntime),WKM(Ntime))
      Allocate (Coff(Nwave,M), Sf(Nwave,MP1))              !, Freqf(Nwave,MP1))
      Allocate (SMean(Nwave), Dfre(Nwave), Period(MP1))    !, Period(Nwave,MP1))

     ! Calculate Frequency fdt=f*Delta_S and Period:
      Do I=1,MP1
         Freq(I) = (I-1)*DELTA_S
         if (I /= 1) then
            Period(i) = 1./Freq(I)
         else
            Period(i) = 2*Ntime !/2.0 ! Miss
         endif
      Enddo
     ! Calculate Power Spectrum over all Nwave:      
      Do iw=1, Nwave
         print*,"Calculating mode ",iw,"/",Nwave
         DAT1(1:Ntime) = imf(1:Ntime,iw)      ! Data of each series
         CALL MEMCOF(DAT1,Ntime,M,PM,COF)     !,WK1,WK2,WKM)
         Coff(iw,1:M) = COF(1:M)              ! For further use
         
         Stb = 0.0
         Do I=1,MP1
         !   print*,I,MP1
         !   Freq(I) = (I-1)*DELTA_S
            S(I)=EVLMEM(FREQ(I),COF,M,PM)     ! Power Spectrum
            Sf(iw,I) = S(I)
         !   Freqf(iw,i) = Freq(I)
            Stb = Stb + S(I)
         Enddo

         Stb = Stb/FLOAT(MP1)          
         DgrFree = (2*N-0.5*M)/FLOAT(M)
         SMean(iw) = Stb; Dfre(iw) = DgrFree
      
         !DO I=1,MP1
         !   IF (I.NE.1) THEN 
         !       PRI=1./FREQ(I)
         !   ELSE
         !       PRI=Miss
         !   ENDIF    
         !   Period(iw,i) = PRI
         !ENDDO
      Enddo
      Call Write_Spec_to_NC
      Call Write_Spec_to_TXT
!------------------------------------------------------------------------------
      CONTAINS
        Include "/home/tanpv/Mylibs/Sub_MEMCOF.f90"
        include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
        include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
!------------------------------------------------------------------------------
      Subroutine Write_Spec_to_TXT
        Integer :: iw, ifre
        Real, Allocatable :: Xwave(:), Xfre(:), Peri(:), tmp

        Allocate (Xwave(1:Nwave),Xfre(1:MP1),Peri(1:MP1))
        Do iw=1, Nwave; Xwave(iw)=iw; Enddo
        Do ifre=1, MP1
        !   IF (Ifre/=1) THEN
              Peri(ifre) = Period(ifre) ! 1./((ifre-1)*DELTA_S)
              Xfre(ifre) = Freq(ifre)   ! (ifre-1)*DELTA_S
        !   Else
        !      Peri(ifre)=999
        !      Xfre(ifre) = 0
        !   Endif
        Enddo
        Fname="Spectral.txt"
        Open (3, file=trim(Fname), Status="Unknown")
        write (3,"(3A8,20I10)") "No","Freq","Peri",(iw, iw=1,Nwave-1)
        Do ifre=MP1, 1, -1 ! 1, MP1
           write (3,"(I6,2x,F8.5,F8.1,20F10.3)") MP1-ifre+1,Xfre(ifre),Peri(ifre),(Sf(iw,ifre), iw=2,Nwave)
        Enddo
        Close (3)
        Deallocate (Xwave,Xfre,Peri)
      End Subroutine Write_Spec_to_TXT
!------------------------------------------------------------------------------
      Subroutine Write_Spec_to_NC
        Integer :: iw, ifre
        Real, Allocatable :: PSD(:,:), Xwave(:), Xfre(:), Peri(:)

        Allocate (PSD(Nwave,MP1), Xwave(1:Nwave),Xfre(1:MP1),Peri(1:MP1))
        Do iw=1, Nwave
           Do ifre=1, MP1
              PSD(iw,MP1-ifre+1)=Sf(iw,ifre)
           Enddo
        Enddo
        Do iw=1, Nwave; Xwave(iw)=iw; Enddo
        Do ifre=1, MP1
           Peri(MP1-ifre+1)=Log(Period(ifre)) ! Log10(1./((ifre-1)*DELTA_S))
           Xfre(ifre) = Freq(ifre)              ! (ifre-1)*DELTA_S
        Enddo
        Fname="Spectral.nc"
        Call wr_ncdf_2d (trim(Fname),Xwave,Peri,Nwave,MP1,PSD,"PSD", &
             "Power Spectrum","x^2", "mode", "Period", Miss)
!        Call wr_ncdf_2d (trim(Fname),Xwave, Xfre,Nwave,MP1, Sf,"PSD", &
!             "Power Spectrum","x^2", "mode", "Frequency", Miss)

        Call wr_ncdf_1d_apd (trim(Fname), SMean, "Smean",Nwave, &
             "Mean Power Spectrum", "x^2", "mode", Miss)
        Call wr_ncdf_1d_apd (trim(Fname),  Dfre, "Dfre", Nwave, & 
             "Degree of Freedom", "-", "mode", Miss)
        Deallocate (PSD,Xwave,Xfre,Peri)
      End Subroutine Write_Spec_to_NC
!------------------------------------------------------------------------------
      Subroutine Read_Data
        Integer :: it, iw, Nt, Nw
        Real, Allocatable :: Xtime(:), Xwave(:)
        Character (Len=150) :: Fname
         Fname="../EEMD/EEMD_rlts.nc"
         Call rd_ncdf_2d (trim(Fname),Xtime,Xwave,Nt,Nw, imf, "imf01")
         Ntime=Nt; Nwave=Nw
      End Subroutine Read_Data
!------------------------------------------------------------------------------

      END