! 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(:) ! For TC: Real, Allocatable :: Coff1(:), Sf1(:) Integer :: I, iw ! ---------------------------------------------------------------------------------- Real, Allocatable :: TC_NWP(:), Ann_Mean(:), Std(:) ! ---------------------------------------------------------------------------------- !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) Call Read_TC_NC 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)) ! For TC: Allocate ( Period(MP1)) Allocate (Coff1(M), Sf1(M1)) ! 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: CALL MEMCOF(TC_NWP,Ntime,M,PM,COF) !,WK1,WK2,WKM) print*,"Calculating done 1" Coff1(1:M) = COF(1:M) ! For further use print*,COF(1:M) Stb1 = 0.0 Do I=1,MP1 print*,FREQ(I) S(I)=EVLMEM(FREQ(I),COF,M,PM) ! Power Spectrum print*,"i,S(I)=",i,S(I) Sf1(I) = S(I) Stb1 = Stb1 + S(I) Enddo Stb1 = Stb1/FLOAT(MP1) DgrFree1 = (2*N-0.5*M)/FLOAT(M) print*,"Calculating done" !SMean(iw) = Stb; Dfre(iw) = DgrFree !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_TC_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_TC_to_NC Integer :: iw, ifre Real, Allocatable :: PSD(:), Xfre(:), Peri(:), xlon(:) Allocate (PSD(MP1), Xfre(1:MP1),Peri(1:MP1), xlon(MP1)) Do ifre=1, MP1 PSD(MP1-ifre+1)=Sf1(ifre) Enddo Do ifre=1, MP1 Peri(MP1-ifre+1)=Log(Period(ifre)) Xfre(ifre) = Freq(ifre) xlon(ifre) = ifre Enddo Fname="Spectral_TC.nc" Call wr_ncdf_1d (Fname, xlon, MP1, PSD, "spec", & "spectral density", "-", "stt", Miss) Call wr_ncdf_1d_apd (trim(Fname), Xfre, "fre", MP1, & "Frequency", "-", "stt", Miss) Call wr_ncdf_1d_apd (trim(Fname), Peri, "peri", MP1, & "period", "-", "stt", Miss) End Subroutine Write_Spec_TC_to_NC !------------------------------------------------------------------------------ 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 !---------------------------------------------------------------------------------------------------- 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" Call rd_ncdf_2d (Fname,Year,Mon,Nyear,Nmon, TC, "jma_ves_tc") !"JMA_NWP_TC_All") 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 !------------------------------------------------------------------------------ END