! Chuong trinh tinh EOF cua 14 ClimIDX cho cac Lag tu 3-12 thang Use NetCDF ! Su dung SL cua cac thang tu 1979 den dau 2010 de Du bao cho cac thang tu 1/1980 den 12/2010 (Training) ! Real, Parameter :: Miss=-99.0 Integer, Parameter :: Yr11=1980, Yr1=1981, Yr2=2010, Nyear=Yr2-Yr1+1, Nmon=12, Ntime=Nyear*Nmon, Ntime_BL=Ntime+Nmon Integer, Parameter :: Yf11=2009, Yf1=2010, Yf2=2020, Nyfc= Yf2-Yf1+1, Ntime_fc=(Nyfc-1)*Nmon, Nlead=6 Integer, Parameter :: Nt0=(Yf1-Yr1)*Nmon+1 Integer, Parameter :: Nlag=6, Nlag0=2 ! Tre 6 thang, bat dau tu thang thu Nlag0 ke tu thang Du bao (VD: Leadtime=0 tai thang 4, thi Nlag0 tu thang 2) Integer, Parameter :: Nlen=4 Character (Len=7), Parameter :: Dr_nm(Nlen)=(/"spei_01","spei_03","spei_06","spei_12"/) ! Nidx=24: !Integer, Parameter :: Nidx=24 !Character (Len=15), Parameter :: Clim_nm(Nidx)=(/"nino12", "nino3", "nino34", "nino4", "oni","tni",& ! "mei", "censo", "ammsst", "NTA_ersst", & ! "soi", "reqsoi", "whwp", "amonus", "repac_slpa", & ! "rindo_slpa", "tsa", "tna","zwnd20", "wpac85", "epac85", & ! "cpac85", "qbo30", "wp"/) ! Cac chi so co do tre thoi gian 1 thang: ! Nino12, nino3, nino34, nino4, pna, best (censo), ammsst, nta_sst, soi, whwp, amonus, tna, qbo, wp, ao, nao, ReqSOI, ! Repac_slpa, Rindo_slpa, mei_v2, zwnd200, wpac850 !Integer, Parameter :: Nidx=22 !Character (Len=15), Parameter :: Clim_nm(Nidx)=(/"nino12", "nino3", "nino34", "nino4", "pna", & !5 ! "censo", "ammsst", "NTA_ersst", & !3 ! "soi", "whwp", "amonus", "tna", "qbo30", "wp", "ao", & !8 ! "nao", "reqsoi", "repac_slpa", & !3 ! "rindo_slpa", "mei", "zwnd20", "wpac85"/) ! , "epac85", & !4 Integer, Parameter :: Nidx=10 Character (Len=15), Parameter :: Clim_nm(Nidx)=(/"nino3", "nino34", "nino4", "pna", & !5 "censo", & "soi", & "reqsoi", "repac_slpa", & !3 "rindo_slpa", "mei"/) ! , "epac85", & !4 ! "cpac85", "qbo50"/) !2 Real, Allocatable :: Clim_Dat(:,:), Y_Clim(:,:) Real, Allocatable :: Clim_Fcs(:,:), Y_Clim_Fcs(:,:) Character (Len=200) :: FnameTC, FnameIDX, Fname, Out_Dir, Predictand !-------------------------------------------------------------------------------------- ! Tinh cac EOF cua cac ClimIDX va tuong quan giua cac EOF voi TC_Char Real*8, Allocatable :: xdat(:,:), ydat(:,:), xtb(:), sx(:) Real*8, Allocatable :: EOFs(:,:), PCs(:,:) Real*8, Allocatable :: Eigenvalues(:), ECumulate(:) Integer :: Nt, Nx ! Kich thuoc thuc cua mang X Integer :: Mcount, Ncount, switch, Standard, PC_Stard, Nmode, K Real*8 :: Factor Real, Allocatable :: Xtb_Clim(:), Sx_Clim(:), PCs_Clim(:,:), Eig_Clim(:) Real, Allocatable :: lamda(:,:), PC1(:,:) Real, Allocatable :: EOF_Clim(:,:), ClimIDX(:,:) ! (Nidx,Ntime+Nmon) Loai bo CPAC850 va WPC850 vi SL ngan Real, Allocatable :: EOF_Fcs(:,:), EOF1(:,:), ClimFcs(:,:) Integer :: imon, iyr, iidx, imod, Nt_Fcs Real, Allocatable :: Corr_all_lag(:,:,:),Prob_all_lag(:,:,:) ! Nmon,Nmode,NTC_Cha !-------------------------------------------------------------------------------------- ! So bien duoc chon boi Stepwise Regression Integer, Allocatable :: Mvar (:,:) ! (NTC_Cha,Nmon) Integer, Allocatable :: IDX_Var(:,:,:) ! (NTC_Cha,Nx,Nmon) Integer, Parameter :: Mvar_X=9 Real, Allocatable :: X_Clim(:,:,:,:), X_Fcs(:,:,:,:) Integer :: NMax, MMax, N, M Real :: PCT, ANS(11) Integer :: NSTEP(5), i,j, NN, KSTEP Integer, Allocatable :: IDX(:),L(:) Real, Allocatable :: X(:),XBAR(:),STD(:),D(:),B(:),RX(:),R(:),T(:) Integer, Allocatable :: NSTEPS(:,:),LS(:,:) Real, Allocatable :: BS(:,:), ANSS(:,:) Real, Allocatable :: XMEAN(:),DX(:), RR(:,:), YORG(:), YHQ(:) !-------------------------------------------------------------------------------------- Real, Allocatable :: DrIDX(:,:,:), Dr_BL(:,:,:,:), Dr_Fcs(:,:,:,:) Real, Allocatable :: Lon(:), Lat(:) Character (Len=20),Allocatable :: Stname(:) Integer :: Nst, ist, Ntime_dr !-------------------------------------------------------------------------------------- Real, Allocatable :: Beta(:,:,:,:) ! Cac he so HQ Real, Allocatable :: Yout(:,:,:,:) ! Y-estimated Real, Allocatable :: Yfcs(:,:,:,:) ! Y-estimated for Independent forecast Real, Allocatable, Dimension(:,:,:) :: ME0, MAE0, RMSE0, RME0, Corr0, Pvalue0 ! Training Real, Allocatable, Dimension(:,:,:) :: ME1, MAE1, RMSE1, RME1, Corr1, Pvalue1 ! Forecast !-------------------------------------------------------------------------------------- Predictand="spei" Out_Dir="/home/tanpv/Progs/2022/Dr_VNU/RES/" FnameTC="/home/tanpv/Progs/2022/Dr_VNU/RES/SPEI_Stn_VN_1980-2020.nc" FnameIDX="/work/projects/G09/SEACLID/ClimIDX_2022/Clim_IDX_2022.nc" print*, "Calling Read_ClimIDX" Call Read_ClimIDX (FnameIDX) Call Read_Dr_IDX (FnameTC) print*, "Calling Multi_lag" Call Multi_lag ! Return: Clim_Dat switch=0; Standard=1; PC_Stard=0; Nmode=6 Nt=Ntime; Nx=Nidx*Nlag !-------------------------------------------------------------------------------------- if (switch==1) then Mcount=Nt; Ncount=Nx else if (switch==0) then Mcount=Nx; Ncount=Nt endif Factor = Sqrt(Real(Nt)) Allocate (xdat(Nt,Nx), ydat(Nt,Nx)) Allocate (Xtb(Nx), Sx(Nx)) ! TB, DLC, ma tran Covar hoac Corre Allocate (Eigenvalues(Mcount)) ! Tri rieng va Vecto rieng (1D) Allocate (EOFs(Nt,Nmode), PCs(Nx,Nmode), ECumulate(Nmode)) Allocate (EOF_Clim(Nt,Nmode)) Allocate (Xtb_Clim(Nx), Sx_Clim(Nx), PCs_Clim(Nx,Nmode), Eig_Clim(Nmode)) !-------------------------------------------------------------------------------------- ! PCA: Do iyr=1,Nt !=Ntime Do iidx=1,Nx ! =Nidx*Nlag xdat(iyr,iidx)=Clim_Dat(iidx,iyr) Enddo Enddo Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, & Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode) print"(A,6(I2,F8.2),F10.2)","Eigen - mode ",(imod, Eigenvalues(imod)/sum(Eigenvalues)*100,imod=1,Nmode),sum(Eigenvalues(1:Nmode))/sum(Eigenvalues)*100. Do imod=1,Nmode Do iyr=1,Nt EOF_Clim(iyr,imod)=EOFs(iyr,imod) Enddo Do iidx=1,Nx PCs_Clim(iidx,imod)=PCs(iidx,imod) Enddo Eig_Clim(imod)=Eigenvalues(imod) Enddo Do iidx=1,Nx Xtb_Clim(iidx)=Xtb(iidx); Sx_Clim(iidx)=Sx(iidx) Enddo ! Wtite to file for further forecast: Fname="RES/pca.nc" Call Write_PCs_Xbar_Std_Eig_for_Fcs (Fname) ! Test: !Allocate (Y_Clim(Nt,Nx)) !Do iyr=1,Nt ! Do iidx=1,Nx ! Y_Clim(iyr,iidx)=Xdat(iyr,iidx)-Xtb_Clim(iidx) ! if (Standard==1) Y_Clim(iyr,iidx)=Y_Clim(iyr,iidx)/Sx_Clim(iidx) !/Factor ! Y(Nt,Nx)=(X(Nt,Nx)-Xtb(Nx))/Sx(Nx)/sqrt(Nt) ! Enddo ! print"(3F8.3)",(Y_Clim(iyr,iidx),Ydat(iyr,iidx),Ydat(iyr,iidx)-Y_Clim(iyr,iidx),iidx=1,Nx); pause !Enddo !print*,Xtb_Clim; pause; !print*,Sx_Clim; pause !print*,"EOF_Clim:" !print"(9F10.2)", Matmul(transpose(EOF_Clim),EOF_Clim); pause !print*,"PCs_Clim:" !print"(9F10.2)", Matmul(transpose(PCs_Clim),PCs_Clim); pause ! Calculating EOF_Fcs: Nt_Fcs=Ntime_fc ! Bo mot nam Allocate (Y_Clim_Fcs(Nt_Fcs,Nx),EOF_Fcs(Nt_Fcs,Nmode),EOF1(Nt_Fcs,Nmode)) Allocate (lamda(Nmode,Nmode), PC1(Nx,Nmode)) lamda = 0.0 ! Set all values to zero Do iyr=1,Nt_Fcs ! Calculating anomalous values of Fcs: Do iidx=1,Nx Y_Clim_Fcs(iyr,iidx)=Clim_Fcs(iidx,iyr)-Xtb_Clim(iidx) if (Standard==1) Y_Clim_Fcs(iyr,iidx)=Y_Clim_Fcs(iyr,iidx)/Sx_Clim(iidx) ! Y(Nt,Nx)=(X(Nt,Nx)-Xtb(Nx))/Sx(Nx) !/sqrt(Nt) Enddo !print"(14F8.3)",Y_Clim_Fcs(iyr,1:Nx); pause Enddo !pause !PC1(1:Nx,1:Nmode)=PCs_Clim(1:Nx,1:Nmode) !print"(6F8.3)",matmul(transpose(PC1),PC1); pause Do j=1,Nmode lamda(j,j) = 1./Eig_Clim(j) ! 1/sqrt(eig) Enddo !EOF1 = Matmul(Matmul(Y_Clim_Fcs,PC1),lamda)/Factor ! F(Nt,K) = Y(Nt,Nx) x A(Nx,K) x lamda(K,K) / sqrt(Nt) EOF_Fcs = Matmul(Matmul(Y_Clim_Fcs,PCs_Clim),lamda)/Factor ! F(Nt,K) = Y(Nt,Nx) x A(Nx,K) x lamda(K,K) / sqrt(Nt) ! Test: print*,"EOF_Fcs:" print"(6F10.2)", Matmul(transpose(EOF_Fcs),EOF_Fcs); pause ! print test: Open (3,file="ANN_Inp.txt",status="unknown") ! ist=1, ilt=3, ilen=3 write (3,*) "Input for ANN 1 hidden layer" write (3,"(2I6)") Nt+Nt_Fcs, Nmode+2 ! Nt+Nt_Fcs = Ca chuoi EOF PT va DL; Nmode+2 = Them Bien PT duoc su dung nhu bien DL Do iyr=1,Nt write(3,"(15F8.3)") Dr_BL(1,3,3,iyr), Dr_BL(1,3,3,iyr), EOF_Clim(iyr,:) Enddo ! print* Do iyr=1,Nt_Fcs write(3,"(15F8.3)") Dr_Fcs(1,3,3,iyr), Dr_Fcs(1,3,3,iyr), EOF_Fcs(iyr,:) Enddo Close (3) !-------------------------------------------------------------------------------------- ! Training and Forecast: Allocate (Beta(Nst,Nlead,Nlen,0:Nmode),Yout(Nst,Nlead,Nlen,Ntime)) Allocate (Yfcs(Nst,Nlead,Nlen,Ntime_fc)) !Call Calc_HQTT ! Verification: Allocate (ME0(Nst,Nlead,Nlen),MAE0(Nst,Nlead,Nlen),Corr0(Nst,Nlead,Nlen),Pvalue0(Nst,Nlead,Nlen)) Allocate (ME1(Nst,Nlead,Nlen),MAE1(Nst,Nlead,Nlen),Corr1(Nst,Nlead,Nlen),Pvalue1(Nst,Nlead,Nlen)) Call Verification ! Write out to GrADS Station files: Call Write_GrADS_Stn_Multivar (ME0,Nst,Nlead,Nlen,Lon,Lat,Stname,"me",trim(Out_Dir)//trim(Predictand)//"_ME_Training") Call Write_GrADS_Stn_Multivar (MAE0,Nst,Nlead,Nlen,Lon,Lat,Stname,"mae",trim(Out_Dir)//trim(Predictand)//"_MAE_Training") Call Write_GrADS_Stn_Multivar (Corr0,Nst,Nlead,Nlen,Lon,Lat,Stname,"cor",trim(Out_Dir)//trim(Predictand)//"_COR_Training") Call Write_GrADS_Stn_Multivar (Pvalue0,Nst,Nlead,Nlen,Lon,Lat,Stname,"pval",trim(Out_Dir)//trim(Predictand)//"_Pval_Training") Call Write_GrADS_Stn_Multivar (ME1,Nst,Nlead,Nlen,Lon,Lat,Stname,"me",trim(Out_Dir)//trim(Predictand)//"_ME_Fcs") Call Write_GrADS_Stn_Multivar (MAE1,Nst,Nlead,Nlen,Lon,Lat,Stname,"mae",trim(Out_Dir)//trim(Predictand)//"_MAE_Fcs") Call Write_GrADS_Stn_Multivar (Corr1,Nst,Nlead,Nlen,Lon,Lat,Stname,"cor",trim(Out_Dir)//trim(Predictand)//"_COR_Fcs") Call Write_GrADS_Stn_Multivar (Pvalue1,Nst,Nlead,Nlen,Lon,Lat,Stname,"pval",trim(Out_Dir)//trim(Predictand)//"_Pval_Fcs") !-------------------------------------------------------------------------------------- CONTAINS !-------------------------------------------------------------------------------------- Subroutine ANN_1 (input_file, output_file, save_file) character (len=*) :: input_file, output_file, save_file logical :: training_mode ! Use for calculating scaling factors Integer :: Nsample, Nepoc Real,allocatable, Dimension(:,:) :: Inputs, Outputs Real,allocatable, Dimension(:) :: min_in, max_in, min_out, max_out Real :: mean_error integer :: Nin, Nhidden, Nout, iin, ihid, iout integer :: i,j,k Real,allocatable, Dimension(:,:) :: weight_ih, weight_ho, moment_ih, moment_ho Real,allocatable, Dimension(:) :: fin, fhid, fout, ftarget, delta_o, delta_h ! Training Real :: learning_rate,momentum, error logical :: debug Nhidden = 2 ! Number of hidden layers learning_rate = 0.5 momentum = 0 training_mode = .true. ! For training Nepoc = 8000 ! Iterative training times debug = .false. !input_file = trim(FnameInp(ifile)) !output_file = trim(FnameOut(ifile)) !save_file = trim(FnameSave(ifile)) print*, "Calculating ",input_file Call controller_init if (training_mode) then ! .T. = Training, .F. = forecast (Define from Namelist) Call controller_train_data endif write(*,*) "Writting output..." Call controller_output_data !(Yhq) Call controller_deinit ! Yout(ist,ilt,ilen,1:Ntime)=Yhq(1:Ntime) Yout(1,3,3,1:Ntime)=Yhq(1:Ntime) ! ! ist=1, ilt=3, ilen=3 End Subroutine ANN_1 !-------------------------------------------------------------------------------------- Subroutine Verification ! Subroutine Calc_Verification (Yorg,Yhq,N, ME, MAE, RMSE, RME, Corr, Pvalue) ! Real, Allocatable, Dimension(:,:,:) :: ME, MAE, RMSE, RME, Corre, Pvalue Integer :: ist,ilt,ilen Real :: ME, MAE, RMSE, RME, Corr, Pvalue Real, Allocatable :: Yorg(:),Yhq(:) ! Training: Allocate (Yorg(Ntime),Yhq(Ntime)) Do ilen=1,Nlen Do ist=1,Nst Do ilt=1,Nlead Yorg(1:Ntime)=Dr_BL(ist,ilt,ilen,1:Ntime) Yhq (1:Ntime)=Yout (ist,ilt,ilen,1:Ntime) Call Calc_Verification (Yorg,Yhq,Ntime, ME, MAE, RMSE, RME, Corr, Pvalue) ME0(ist,ilt,ilen)=ME MAE0(ist,ilt,ilen)=MAE Corr0(ist,ilt,ilen)=Corr; Pvalue0(ist,ilt,ilen)=Pvalue Enddo Enddo Enddo Deallocate (Yorg,Yhq) ! Forecast: Allocate (Yorg(Ntime_fc),Yhq(Ntime_fc)) Do ilen=1,Nlen Do ist=1,Nst Do ilt=1,Nlead Yorg(1:Ntime_fc)=Dr_Fcs(ist,ilt,ilen,1:Ntime_fc) Yhq (1:Ntime_fc)=Yfcs (ist,ilt,ilen,1:Ntime_fc) Call Calc_Verification (Yorg,Yhq,Ntime_fc, ME, MAE, RMSE, RME, Corr, Pvalue) ME1(ist,ilt,ilen)=ME MAE1(ist,ilt,ilen)=MAE Corr1(ist,ilt,ilen)=Corr; Pvalue1(ist,ilt,ilen)=Pvalue Enddo Enddo Enddo End Subroutine Verification !-------------------------------------------------------------------------------------- Subroutine Calc_HQTT Real, Allocatable :: X(:),Beta_tmp(:),Yhq(:) ! X(N*(Mvars+1)) Real, Allocatable :: X1(:), Y1(:), Rxy(:) Real :: S,Q,U,RR,F, Cor ! Dr_BL(ist,ilt,ilen,it1), EOF_Clim(iyr,imod)=EOFs(iyr,imod), Mvars=Nmode+1 ! Test correlations: Allocate (X1(Ntime),Y1(Ntime),Rxy(Nmode)) Do ilen=1,Nlen Do ist=1,Nst Do ilt=1,Nlead X1(1:Ntime)=Dr_BL(ist,ilt,ilen,1:Ntime) Do imod=1,Nmode Y1(1:Ntime)=EOF_Clim(1:Ntime,imod) Call Correlation (X1,Y1,Ntime, Cor, Pval, Miss); Rxy(imod)=Cor Enddo ! print"(A,3I4,10F10.3)","Cor ilen,ist,ilt: ",ilen,ist,ilt,Rxy(:) Enddo Enddo Enddo ! PAUSE ! Training: Allocate (X(Ntime+Ntime*Nmode),Beta_tmp(Nmode+1),Yhq(Ntime)) it1=Ntime Do imod=1,Nmode Do it=1,Ntime it1=it1+1; X(it1)=EOF_Clim(it,imod) Enddo Enddo Do ilen=1,Nlen Do ist=1,Nst Do ilt=1,Nlead Do it=1,Ntime X(it)=Dr_BL(ist,ilt,ilen,it) ! Ntime gia tri dau la Y Enddo print*,"Calling REGRE for Len, Stn, LT: ",ilen,ist,ilt Call REGRE(X,Ntime,Nmode+1,Beta_tmp, S,Q,U,RR,F) !print*,"Beta: ",Beta_tmp;pause Do it=1,Ntime Yhq(it)=0 Do j=2,Nmode+1 Yhq(it)=Yhq(it)+Beta_tmp(j)*EOF_Clim(it,j-1) Enddo Yhq(it)=Yhq(it)+Beta_tmp(1) ! +A0 print"(3F10.3)",Dr_BL(ist,ilt,ilen,it),Yhq(it),Yhq(it)-Dr_BL(ist,ilt,ilen,it) Enddo Yout(ist,ilt,ilen,1:Ntime)=Yhq(1:Ntime) ! Estimated Beta(ist,ilt,ilen,0:Nmode)=Beta_tmp(1:Nmode+1) Enddo Enddo Enddo Deallocate(X,Beta_tmp,Yhq) !pause Allocate (Yhq(Ntime_fc)) ! Forecast: Do ilen=1,Nlen Do ist=1,Nst Do ilt=1,Nlead Do it=1,Ntime_fc Yhq(it)=0 Do j=1,Nmode Yhq(it)=Yhq(it)+Beta(ist,ilt,ilen,j)*EOF_Fcs(it,j) Enddo Yhq(it)=Yhq(it)+Beta(ist,ilt,ilen,0) ! +A0 !print"(3F10.3)",Dr_Fcs(ist,ilt,ilen,it),Yhq(it),Yhq(it)-Dr_Fcs(ist,ilt,ilen,it) Enddo Yfcs(ist,ilt,ilen,1:Ntime_fc)=Yhq(1:Ntime_fc) ! Estimated Enddo Enddo Enddo End Subroutine Calc_HQTT !-------------------------------------------------------------------------------------- Subroutine Write_PCs_Xbar_Std_Eig_for_Fcs (Fname) Character (Len=*) :: Fname Integer :: i,j Real, Allocatable :: xvar(:), xmod(:), X1d(:) ! Xtb_Clim(iidx)=Xtb(iidx); Sx_Clim(iidx)=Sx(iidx) Allocate (xvar(Nx),xmod(Nmode),X1d(Nx)) Do i=1,Nx; xvar(i)=i; Enddo Do i=1,Nmode; xmod(i)=i; Enddo Call wr_ncdf_2d (Fname,xvar,xmod,Nx,Nmode,PCs_Clim,"pc","PCs of ClimIDX","-","var","mode",Miss) Call wr_ncdf_1d_apd (Fname,Xtb_Clim,"mean",Nx,"Mean values of Variables","-","var",Miss) Call wr_ncdf_1d_apd (Fname,Sx_Clim,"std",Nx,"Standard deviation of Variables","-","var",Miss) X1d=Miss; X1d(1:Nmode)=Eig_Clim(1:Nmode) Call wr_ncdf_1d_apd (Fname,X1d,"eig",Nx,"Eigen values","-","var",Miss) Deallocate (xvar,xmod,X1d) End Subroutine Write_PCs_Xbar_Std_Eig_for_Fcs !-------------------------------------------------------------------------------------- Subroutine Read_ClimIDX (Fname) Character (Len=*) :: Fname Real, Allocatable :: Year(:),Month(:), tmp(:,:) Integer :: Nyr, Nmn, iidx, iyr, i,j, k, imon !Yr11=Yr1-1 ! Lui lai 1 nam Allocate (ClimIDX(Nidx,Ntime_BL), ClimFcs(Nidx,Ntime_fc+Nmon)) Do iidx=1,Nidx Call rd_ncdf_2d (Fname,Year,Month,Nyr,Nmn,tmp, trim(Clim_nm(iidx))) print*,"Reading ",iidx," ",trim(Clim_nm(iidx)) j=0 Do iyr=Yr11,Yr2 Do i=1,Nyr if (Year(i)==iyr) then Do imon=1,Nmon j=j+1; ClimIDX(iidx,j)=tmp(i,imon) !print "(F5.0,12F8.3)",Year(i) !,ClimIDX(iidx,j) print*,j Enddo endif Enddo Enddo !print "(12F8.3)",ClimIDX(iidx,:); pause ! For Independent Forecast: j=0 Do iyr=Yf11,Yf2 Do i=1,Nyr if (Year(i)==iyr) then Do imon=1,Nmon j=j+1; ClimFcs(iidx,j)=tmp(i,imon) print "(F5.0,12F8.3)",Year(i) !,ClimFcs(iidx,j,:) Enddo endif Enddo Enddo !print*,trim(Clim_nm(iidx));pause Deallocate (Year,Month,tmp) !print*,Clim_nm(iidx);pause Enddo End Subroutine Read_ClimIDX !-------------------------------------------------------------------------------------- Subroutine Multi_lag Integer :: it, imon, iidx, ilag, imon1, it1, it2, k, Nlg0 Integer :: mn,yrs Real, Allocatable :: X(:), Y(:) Nlg0=Nmon-Nlag ! Lt=0 tai 01/1980, ClimIDX tinh tu Nmon-Lag0+1 = 11/1979 Allocate (Clim_Dat(Nidx*Nlag,Ntime),Clim_Fcs(Nidx*Nlag,Ntime_fc)) k=0 Do iidx=1, Nidx ! imon=1, Nmon Do ilag=1,Nlag k=k+1; it1=0 Do it=Nlg0+ilag-1,Nlg0+ilag-2+Ntime it1=it1+1; Clim_Dat(k,it1)=ClimIDX(iidx,it) !print*,k,it1,it !mn=mod(it,Nmon); if (mn==0) mn=Nmon !yrs=it/Nmon; if (yrs*Nmon/=it) yrs=yrs+1; yrs=yrs+Yr11-1 !print"(A,6I6,F10.4)","it,yrs,mn,iidx,ilag,it1=",it,yrs,mn,iidx,ilag,it1,Clim_Dat(k,it1) Enddo !pause Enddo Enddo ! For Independent Forecast: k=0 Do iidx=1, Nidx ! imon=1, Nmon Do ilag=1,Nlag k=k+1; it1=0 Do it=Nlg0+ilag-1,Nlg0+ilag-2+Ntime_fc it1=it1+1; Clim_Fcs(k,it1)=ClimFcs(iidx,it) !mn=mod(it,Nmon); if (mn==0) mn=Nmon !yrs=it/Nmon; if (yrs*Nmon/=it) yrs=yrs+1; yrs=yrs+Yf11-1 !print"(A,6I6,F10.4)","it,yrs,mn,iidx,ilag,it1=",it,yrs,mn,iidx,ilag,it1,Clim_Fcs(k,it1) Enddo !pause Enddo Enddo ! print test: !Do iidx=1,Nidx ! it1=Nlag*(iidx-1)+1; it2=it1+Nlag-1 ! Do it=1,Ntime_fc ! print"(30F6.2)",Clim_Fcs(it1:it2,it) !,Clim_Dat(23:28,it),Clim_Dat(45:50,it) ! Enddo ! print*,iidx," ", Clim_nm(iidx); pause !Enddo ! DrIDX for training: !Allocate (Dr_BL(Nst,Ntime,Nlead,Nlen),Dr_Fcs(Nst,Ntime_fc,Nlead,Nlen)) Allocate(Dr_BL(Nst,Nlead,Nlen,Ntime),Dr_Fcs(Nst,Nlead,Nlen,Ntime_fc)) Do ist=1,Nst ! print*,"Multi lad for station ",ist Do ilt=1,Nlead ! For training period Do ilen=1,Nlen it1=0 Do it=ilt,ilt+Ntime-1 it1=it1+1 !Do ilen=1,Nlen !Dr_BL(ist,it1,ilt,ilen)=DrIDX(ist,it,ilen) Dr_BL(ist,ilt,ilen,it1)=DrIDX(ist,it,ilen) !mn=mod(it,Nmon); if (mn==0) mn=Nmon !yrs=it/Nmon; if (yrs*Nmon/=it) yrs=yrs+1; yrs=yrs+Yr1-1 !print"(A,7I6,F10.4)","it,yrs,mn,ist,ilt,ilen,it1=",it,yrs,mn,ist,ilt,ilen,it1,Dr_BL(ist,ilt,ilen,it1) Enddo !pause Enddo Enddo Enddo !pause Do ist=1,Nst Do ilt=1,Nlead ! For forecating period Do ilen=1,Nlen it1=0 Do it=Nt0+ilt-1,Nt0+ilt+Ntime_fc-2 it1=it1+1 !Do ilen=1,Nlen !Dr_Fcs(ist,it1,ilt,ilen)=DrIDX(ist,it,ilen) Dr_Fcs(ist,ilt,ilen,it1)=DrIDX(ist,it,ilen) !mn=mod(it,Nmon); if (mn==0) mn=Nmon !yrs=it/Nmon; if (yrs*Nmon/=it) yrs=yrs+1; yrs=yrs+Yr1-1 !print"(A,7I6,F10.4)","it,yrs,mn,ist,ilt,ilen,it1=",it,yrs,mn,ist,ilt,ilen,it1,Dr_Fcs(ist,ilt,ilen,it1) Enddo !pause Enddo Enddo Enddo !pause End Subroutine Multi_lag !-------------------------------------------------------------------------------------- Subroutine Read_Dr_IDX (Fname) Character (Len=*) :: Fname Real, Allocatable :: Stn(:),xmon(:), tmp(:,:) Integer :: Nstn, Nmn, icha, iyr, i,j, it1, it2 ! DrIDX, Dr_BL(:,:,:), Dr_Fcs(:,:,:) Call rd_ncdf_2d (Fname,Stn,xmon,Nstn,Nmn,tmp,"spei_01") Nst = Nstn ! Num of Stn it1=Nmon+1; it2=Nmn ! Bo nam dau tien Ntime_dr = it2-it1+1 Allocate (DrIDX(Nst,Ntime_dr,Nlen)) print*,"Nst,Ntime_dr,Nlen=",Nst,Ntime_dr,Nlen;pause DrIDX(:,1:Ntime_dr,1) = tmp(:,it1:it2); Deallocate (Stn,xmon,tmp) Call rd_ncdf_2d (Fname,Stn,xmon,Nstn,Nmn,tmp,"spei_03") DrIDX(:,1:Ntime_dr,2) = tmp(:,it1:it2); Deallocate (Stn,xmon,tmp) Call rd_ncdf_2d (Fname,Stn,xmon,Nstn,Nmn,tmp,"spei_06") DrIDX(:,1:Ntime_dr,3) = tmp(:,it1:it2); Deallocate (Stn,xmon,tmp) Call rd_ncdf_2d (Fname,Stn,xmon,Nstn,Nmn,tmp,"spei_12") DrIDX(:,1:Ntime_dr,4) = tmp(:,it1:it2); Deallocate (Stn,xmon,tmp) Call rd_ncdf_1d (Fname,Stn,Nstn,Lon,"lon"); Deallocate (Stn) Call rd_ncdf_1d (Fname,Stn,Nstn,Lat,"lat"); Deallocate (Stn) Call rd_ncdf_1d_Char (Fname,Stn,Nstn,Stname,"stname"); Deallocate (Stn) End Subroutine Read_Dr_IDX !-------------------------------------------------------------------------------------- Subroutine Calc_Verification (Yorg,Yhq,N, ME, MAE, RMSE, RME, Corr, Pvalue) ! Input: Yorg(1:N) is an array stored the origin values of Y ! Xorg(1:N,1:M) is an array stored the origin values of X ! Beta(0:M) is an array stored the regression coeficients B0, B1,..., BM ! Output: ME is Mean Error ! MAE is Mean Absolute Error ! RMSE is Root Mean Square Error ! RME is Relative Mean Error ! Corr, Pvalue are Correlation Coefficient and P-value (t-test) Integer :: N, i, j, k, Nc, Nc0 Real :: Yorg(N), Yhq(N), Error Real :: ME, MAE, RMSE, RME, Corr, Pvalue ME=0; MAE=0; RMSE=0; RME=0; Nc=0; Nc0=0 Do i=1, N if (Yhq(i)/=Miss.AND.Yorg(i)/=Miss) then Nc=Nc+1 Error= Yhq(i)-Yorg(i) ME = ME + Error MAE = MAE + ABS(Error) RMSE = RMSE + Error*Error if (Yorg(i)/=0) then; Nc0=Nc0+1; RME = RME + ABS(Error/Yorg(i)); endif endif Enddo ME = ME/Real(Nc); MAE = MAE/Real(Nc); RME = RME/Real(Nc0)*100.; RMSE = Sqrt(RMSE/Real(Nc)) Call Correlation (Yorg, Yhq, N, Corr, Pvalue, Miss) End Subroutine Calc_Verification !-------------------------------------------------------------------------------------- Subroutine Write_GrADS_Stn_Multivar (X,Nst,Nvar,Nt,xlon,xlat,Stname,varn,Grads_File) ! Inputs: X(Nst,Nvar,Nt) - 3D array ! varn - VarName ! Grads_File - Output Filename, including directory Integer :: imon, ivar, ist, ic, i, iidx Integer :: Nst,Nvar,Nt ! Station, Var, Time Character (Len=8) :: STID Character *(*) :: varn,Grads_File Real :: X(Nst,Nvar,Nt) ! X(Nst,Nmode,Nlen,ndr_idx) !Nst,Nmon,NLen,ndr_idx) ! (Nst,2,Nvar) Real :: xlon(Nst), xlat(Nst) Character (Len=20) :: Stname(Nst) Real :: RLAT, RLON, TIM, tmp, Var Integer :: NLEV,NFLAG !---------------------------------------------------------------------------------- Open (10,file=trim(Grads_File)//'.dat',FORM='UNFORMATTED',RECORDTYPE='STREAM') TIM = 0.0 ! Co the bien doi tu -0.5 .. +0.5; thuong dat bang 0. NFLAG = 1 ! Khong doi Do i = 1,Nt !Nlen ! Nvar chi so, moi chi so duoc xem nhu 1 timeslice NLEV = 1 ! =0 o ban ghi cuoi cung cua 1 timeslice; =1 cho tat ca cac timeslice khac ic = 0 Do ist=1, Nst ic = ic+1 STID = Stname(ist)(1:8) ! Chi lay 8 ky tu dau tien RLAT = xlat(ist); RLON = xlon(ist) WRITE (10) STID,RLAT,RLON,TIM,NLEV,NFLAG Do ivar=1,Nvar !Nmode Var = X(ist,ivar,i) WRITE (10) Var ! Gia tri cua bien Enddo Enddo NLEV = 0 WRITE (10) STID,RLAT,RLON,TIM,NLEV,NFLAG Enddo ! ivar Close(10) !---------------------------------------------------------------------! ! Write control file Station data ! !---------------------------------------------------------------------! open (2001,file=trim(Grads_File)//'.ctl') write (2001,'(2A)') 'dset ','^'//trim(Grads_File)//'.dat' write (2001,'(A)') 'dtype station' write (2001,'(2A)') 'stnmap ','^'//trim(Grads_File)//'.dat'//'.map' write (2001,'(A,F8.1)') 'undef ',Miss write (2001,'(A)') 'title Factor Analysis' write (2001,'(A,I3,1x,A)') 'tdef ',Nt,' linear jan1980 1mo' write (2001,'(A,I3)') 'vars ',Nvar Do ivar=1,Nvar write (2001,'(A,I2.2,A)') trim(varn),ivar,' 0 99 Drough Chars' Enddo write (2001,'(A)') 'endvars' close (2001) End Subroutine Write_GrADS_Stn_Multivar !-------------------------------------------------------------------------------------- Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Sub_PCA_2019.f90" Include "/home/tanpv/Mylibs/Calc_Correlation_Coefficient.f90" Include "Sub_neural.f90" !Include "/home/tanpv/Mylibs/Sub_HQTT.f90" !Include "/home/tanpv/Mylibs/Sub_HQTB.f90" !-------------------------------------------------------------------------------------- END !-------------------------------------------------------------------------------------- !Include "/home/tanpv/Mylibs/Sub_HQTB.f90" Include "/home/tanpv/Mylibs/Sub_HQTT.f90" Include "/home/tanpv/Mylibs/Student_Prob.f90" !--------------------------------------------------------------------------------------