PROGRAM runEEMD Use NetCDF IMPLICIT NONE INCLUDE 'eemd.common' !INTEGER, PARAMETER :: MAXSIZE=1000000 !INTEGER, PARAMETER :: IMFS=12 REAL, PARAMETER :: STDEV=0.2 INTEGER, PARAMETER :: EMEMBERS=100 Real, Parameter :: Miss=-99.0 !! INTEGER, PARAMETER :: IMFS=12 INTEGER :: i,j,TimeSeriesLength REAL, DIMENSION(MAXSIZE) :: indata INTEGER:: idum REAL :: eemd_result(MAXSIZE,IMFS+2) Real, Allocatable :: R_CH(:), EMMD(:,:), Xtime(:), Xwave(:) Logical :: status Character (Len=150) :: Fname idum=3571 Call Read_Rain_NC open (4, file='tempTS.dat',FORM='UNFORMATTED',ACCESS='DIRECT', RECL=1) OPEN (15,FILE='tempEEMD.dat',FORM='UNFORMATTED',ACCESS='DIRECT', RECL=1) ! OPEN(UNIT = 15, FILE = 'newfort.txt', STATUS = 'UNKNOWN') open (19, file='tempEEMD.ctl') do i = 1, MAXSIZE read(4, rec=i, ERR=10) indata(i) print*,i," ",indata(i) enddo 10 continue !! TimeSeriesLength = i-1 !! CALL eemd(TimeSeriesLength,indata,STDEV,EMEMBERS,IMFS,idum,eemd_result) Call eemd(TimeSeriesLength, R_CH, STDEV,EMEMBERS,IMFS,idum,eemd_result) DO i=1,TimeSeriesLength DO j=1,IMFS+2 WRITE(15,REC=i+(j-1)*TimeSeriesLength) (eemd_result(i,j)) ENDDO print "(I5, 30F6.3)", i, (eemd_result(i,j), j=1,IMFS+2) ENDDO 500 FORMAT(14F12.8) WRITE (19,*) 'DSET ^tempEEMD.dat' WRITE (19,*) 'TITLE Sample Model Data' WRITE (19,*) 'UNDEF -9.99E+30 ' WRITE (19,*) 'XDEF 1 linear 0 1' WRITE (19,*) 'YDEF 1 linear 0 1' WRITE (19,*) 'EDEF 14 NAMES TS 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15' WRITE (19,*) 'ZDEF 1 linear 1 1' WRITE (19,*) 'TDEF ',TimeSeriesLength,' LINEAR 1JAN1980 1dys' WRITE (19,*) 'VARS 1 ' WRITE (19,*) ' imf 0 99 100 Data Points' WRITE (19,*) 'ENDVARS' WRITE (19,*) ' ' !----------------------------------------------------------------------- Fname="EEMD_rlts.nc" 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 Call wr_ncdf_2d (trim(Fname),Xtime,Xwave,TimeSeriesLength,IMFS+2,EMMD, & "imf", "Wave", "-", "day", "No", Miss) CONTAINS !---------------------------------------------------------------------------------------------------- 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" !---------------------------------------------------------------------------------------------------- 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 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 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 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 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 runEEMD