PROGRAM Main_EEMD Use NetCDF IMPLICIT NONE INTEGER, PARAMETER :: MAXSIZE=1000000 INTEGER, PARAMETER :: IMFS=5, Nmonth=12 Integer, Parameter :: YrBL1=1961, YrBL2=2019, Nidx=16 Character (Len=10), Parameter :: IDXnm(Nidx)=(/"ONI","MEI","REQSOI","SOI","RINDO","Nino34","CPAC850","BEST", & "WPC850","PDO","PNA","DMI","WHWP","TNI","WP","QBO"/) REAL, PARAMETER :: STDEV=0.2 INTEGER, PARAMETER :: EMEMBERS=100 Real, Parameter :: Miss=-99.0 INTEGER :: i,j,TimeSeriesLength REAL, DIMENSION(MAXSIZE) :: indata INTEGER:: idum REAL :: eemd_result(MAXSIZE,IMFS+2) Real, Allocatable :: R_CH(:), R_CH_mean(:), EMMD(:,:), Xtime(:), Xwave(:) ! CH Rain Real, Allocatable :: TC_NWP(:), Ann_Mean(:), Std(:), ClimIDX(:,:) ! For tesing: Real, Allocatable :: xdump(:) ! Timeseries created by sum of different sin functions ! Logical :: status Character (Len=150) :: Fname Integer :: iidx !---------------------------------------------------------------------------------------------------- idum=3571 Fname="EEMD_rlts.nc" !----------------------------------------------------------------------- !Call Read_Rain_NC !Call eemd(TimeSeriesLength, R_CH, STDEV,EMEMBERS,IMFS,idum,eemd_result) !Call Write_EEMD_to_NC(Fname,1) !----------------------------------------------------------------------- Call Read_TC_NC Call eemd(TimeSeriesLength, TC_NWP, STDEV,EMEMBERS,IMFS,idum,eemd_result) Call Write_EEMD_to_NC(Fname,1) !----------------------------------------------------------------------- !Call Read_ClimIDX !Do iidx=1, Nidx ! Call eemd(TimeSeriesLength,ClimIDX(:,iidx),STDEV,EMEMBERS,IMFS,idum,eemd_result) ! Call Write_EEMD_to_NC(Fname,iidx) !Enddo !----------------------------------------------------------------------- ! Test with dummy data: !Call Dummy_Data CONTAINS !---------------------------------------------------------------------------------------------------- Subroutine Dummy_Data Integer :: i,j,k, N, M Real, Allocatable, Dimension (:) :: x1, x2, x3, x4, x5 Real, parameter :: pi=atan(1.) N=500; M=5 Allocate (xdump(N),x1(N), x2(N), x3(N), x4(N), x5(N)) Do i=1,N x1(i)=3.*sin(i*pi/4) x2(i)=0.5*cos(2*i*pi/3) x3(i)=5.*sin(2*pi/5)+0.5*cos(i*pi/7) x4(i)=sin(i*3.)+cos(i*pi/10) x5(i)=2*cos(i*pi/15)+sin(i*pi/2) xdump(i)=x1(i)+x2(i)+x3(i)+x4(i)+x5(i) Enddo Call eemd(N, xdump, STDEV,EMEMBERS,M,idum,eemd_result) ! Write to file Allocate (EMMD(N,M+2),Xtime(N),Xwave(M+2)) Do i=1,N; Xtime(i)=i; Enddo Do j=1,M+2; Xwave(j)=j; Enddo Do i=1,N Do j=1,M+2; EMMD(i,j)=eemd_result(i,j); Enddo Enddo Fname="EEMD_test_with_dummy_data.nc" Call wr_ncdf_2d (trim(Fname),Xtime,Xwave,N,M+2,EMMD, & "imf01", "Wave", "-", "freq", "No", Miss) Call wr_ncdf_1d_apd (trim(Fname),x1,"x1",N,"3.*sin(i*pi/4)","-","freq", Miss) Call wr_ncdf_1d_apd (trim(Fname),x2,"x2",N,"0.5*cos(2*i*pi/3)","-","freq", Miss) Call wr_ncdf_1d_apd (trim(Fname),x3,"x3",N,"5.*sin(2*pi/5)+0.5*cos(i*pi/7)","-","freq", Miss) Call wr_ncdf_1d_apd (trim(Fname),x4,"x4",N,"sin(i*3.)+cos(i*pi/10)","-","freq", Miss) Call wr_ncdf_1d_apd (trim(Fname),x5,"x5",N,"2*cos(i*pi/15)+sin(i*pi/2)","-","freq", Miss) End Subroutine Dummy_Data !---------------------------------------------------------------------------------------------------- Subroutine Write_EEMD_to_NC(Fname,ivar) Character (Len=*) :: Fname Integer :: ivar Character (Len=5) :: varn if (.not.Allocated(EMMD)) Allocate (EMMD(TimeSeriesLength,IMFS+2)) if (.not.Allocated(Xtime)) Allocate (Xtime(TimeSeriesLength)) if (.not.Allocated(Xwave)) Allocate (Xwave(IMFS+2)) ! Allocate (EMMD(TimeSeriesLength,IMFS+2), Xtime(TimeSeriesLength), Xwave(IMFS+2)) Do i=1, TimeSeriesLength; Xtime(i)=i; Enddo Do j=1, IMFS+2; Xwave(j)=j; Enddo Do i=1, TimeSeriesLength Do j=1, IMFS+2 EMMD(i,j)=eemd_result(i,j) Enddo Enddo if (ivar==1) then Call wr_ncdf_2d (trim(Fname),Xtime,Xwave,TimeSeriesLength,IMFS+2,EMMD, & "imf01", "Wave", "-", "day_or_mon_or_Freq", "No", Miss) else write (varn,"(A,I2.2)") "imf",ivar Call wr_ncdf_2d_apd (trim(Fname), EMMD, trim(varn),TimeSeriesLength,IMFS+2, & "Wave", "-", "day_or_mon_or_Freq", "No", Miss) endif End Subroutine Write_EEMD_to_NC !---------------------------------------------------------------------------------------------------- include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" include "/home/tanpv/Mylibs/DayOfMon_NoLeap.f90" include "SUB/eemd.f90" ! include "SUB/gasdev.f90" include "SUB/min_max.f90" include "SUB/natural_spline.f90" include "SUB/standev.f90" !---------------------------------------------------------------------------------------------------- Subroutine Read_ClimIDX Integer :: Nyear, Nmon, iyr, imon, it, iidx Real, Allocatable :: IDX(:,:), Year(:), Mon(:) Character (Len=150) :: Fname Fname="/home/tanpv/OBS/KC09/Clim_IDX.nc" Do iidx=1, Nidx Call rd_ncdf_2d (Fname,Year,Mon,Nyear,Nmon, IDX, trim(IDXnm(iidx))) if (iidx==1) Allocate(ClimIDX(Nyear*Nmon,Nidx)) it=0 Do iyr=1, Nyear Do imon=1, Nmon it=it+1; ClimIDX(it,iidx)=IDX(iyr,imon) Enddo Enddo Deallocate (IDX,Year,Mon) Enddo TimeSeriesLength=Nyear*Nmon End Subroutine Read_ClimIDX !---------------------------------------------------------------------------------------------------- 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") TimeSeriesLength=Nyear*Nmon Allocate (TC_NWP(TimeSeriesLength), Ann_Mean(Nmonth), Std(Nmonth)) ! 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*,Ann_Mean !print*,Std ! Moving summation of Anomaly: !Do it=1, TimeSeriesLength-Mave+1 ! tmp=0 ! Do imon=it, it+Mave-1 ! tmp=tmp+TC_NWP(imon) ! Enddo ! TC_NWP(it)=tmp !Enddo !TimeSeriesLength=TimeSeriesLength-Mave+1 End Subroutine Read_TC_NC !---------------------------------------------------------------------------------------------------- Subroutine Read_Rain_NC Integer, Parameter :: NdayOfYrs=365 ! Real, Parameter :: Miss=-99.0 Integer :: Nst, Nst1, Nyear, Nyrs, Nmon, Nday Integer :: ist, ist1, iyr, iy1, imon, iday, id, id1, id2 Character (Len=150) :: Fname Real, Allocatable :: R(:,:,:,:), Stn(:), Year(:), Mon(:), Day(:), Reg(:) Real, Allocatable :: CH_All(:,:,:), CH(:,:) Real, Allocatable :: Dmean(:,:), X(:) Real :: tmp, Thre Fname="/home/tanpv/OBS_Stn/Daily_NC/R_Daily.nc" Call rd_ncdf_4d (Fname,Stn,Year,Mon,Day,Nst,Nyear,Nmon,Nday,R,"R") Deallocate (Stn) Call rd_ncdf_1d (Fname, Stn, Nst, Reg, "reg") ist1=0 Do ist=1, Nst if (Reg(ist)==6) ist1=ist1+1 Enddo Nst1=ist1; iy1=0 Do iyr=1, Nyear ! Xac dinh nam bat dau cua chuoi SL if (Year(iyr)==1980) then iy1=iyr; exit endif Enddo Allocate (CH_All(Nst1,Nyear,NdayOfYrs),Dmean(Nst1,NdayOfYrs)) ist1=0 ! Luu SL cua CH vao mang CH_All theo cac ngay trong nam Do ist=1, Nst if (Reg(ist)==6) then ist1=ist1+1; Nyrs=0 Do iyr=iy1, Nyear id=0; Nyrs=Nyrs+1 ! Nyrs la so nam cua chuoi se duoc su dung (khac voi Nyear) Do imon=1, Nmon Nday=DayOfMon_NoLeap(imon) Do iday=1, Nday id=id+1 CH_All(ist1,Nyrs,id)=R(ist,iyr,imon,iday) Enddo Enddo Enddo endif Enddo Dmean=Miss; Thre=0.80 ! Tinh trung binh ngay nhieu nam cua tung tram Do ist=1, Nst1 Do id=1, NdayOfYrs iday=0; tmp=0 Do iyr=1, Nyrs if (CH_All(ist,iyr,id) /= Miss) then iday=iday+1; tmp=tmp+CH_All(ist,iyr,id) endif Enddo if (real(iday)/real(Nyrs)>=Thre) Dmean(ist,id)=tmp/real(iday) Enddo Enddo Do ist=1, Nst1 ! Fill up missing data by climatological daily mean Do id=1, NdayOfYrs Do iyr=1, Nyrs if (CH_All(ist,iyr,id)<0.AND.Dmean(ist,id)>0) CH_All(ist,iyr,id)=Dmean(ist,id) Enddo Enddo Enddo Allocate (CH(Nst1,Nyrs*NdayOfYrs), X(Nyrs*NdayOfYrs), R_CH(Nyrs*NdayOfYrs)) Do ist=1, Nst1 ! Daily timeseries data at each station iday=0 Do iyr=1, Nyrs Do id=1, NdayOfYrs iday=iday+1; CH(ist,iday)=CH_All(ist,iyr,id) Enddo Enddo Enddo ist1=0; R_CH=0 Do ist=1, Nst1 ! Areal summation of daily rainfall X(:)=CH(ist,:) if (.Not.Empty_Ar1(X,Nyrs*NdayOfYrs)) then ist1=ist1+1 Do iday=1, Nyrs*NdayOfYrs R_CH(iday)=R_CH(iday)+CH(ist,iday) Enddo endif Enddo Do iday=1, Nyrs*NdayOfYrs ! Areal mean of daily rainfall R_CH(iday)=R_CH(iday)/Real(ist1) Enddo TimeSeriesLength = Nyrs*NdayOfYrs Allocate (R_CH_mean(NdayOfYrs)) ! Calculate Anomaly: Do iday=1, NdayOfYrs id1=0; R_CH_mean(iday)=0 Do id=1, TimeSeriesLength id2=mod(id,NdayOfYrs); if (id2==0) id2=NdayOfYrs if (id2==iday) then id1=id1+1; R_CH_mean(iday)=R_CH_mean(iday)+R_CH(id) endif Enddo R_CH_mean(iday)=R_CH_mean(iday)/Real(id1) Enddo Do id=1, TimeSeriesLength ! id2=mod(id,NdayOfYrs); if (id2==0) id2=NdayOfYrs ! R_CH(id)=R_CH(id)-R_CH_mean(id2) Enddo End Subroutine Read_Rain_NC !------------------------------------------------------------------- Logical Function Empty_Ar1(X,N) Integer :: N, i Real :: X(N) Logical :: Emp Emp=.False. Do i=1, N if (X(i)==Miss) then Emp=.True.; exit endif Enddo Empty_Ar1=Emp End Function Empty_Ar1 !---------------------------------------------------------------------------------------------------- END PROGRAM Main_EEMD include "SUB/gasdev.f90"