! Chuong trinh tinh EOF cua 14 ClimIDX cho cac Lag tu 3-12 thang Use NetCDF Use neural type(network) :: Net Real, Parameter :: Miss=-99.0 Integer, Parameter :: Yr1=1981, Yr2=2020, Nyear=Yr2-Yr1+1, Nmon=12 Integer, Parameter :: Ntime=Nyear*Nmon, Ntime1=Ntime+Nmon, Nlag=2 ! Nlag is used as Nlead Integer, Parameter :: Nlen=2 Integer, Parameter :: Len_SPI(Nlen) = (/6,12/) Integer, Parameter :: Lag(Nlag) = (/4,5/) ! Lag-2 is used as leadtime ; 2 is num of months after predictors are available Integer, Parameter :: Nidx=16 Character (Len=15), Parameter :: Clim_nm(Nidx)=(/"nino3", "nino34", "nino4", "pna", & !5 "censo", "oni", "soi", & "reqsoi", "repac_slpa", & !3 "rindo_slpa", "mei", "dmi", "pdo", & !4 "wpac85","cpac85", "epac85"/) !2 Real, Allocatable :: ClimIDX(:,:) ! (Nidx,Ntime Real, Allocatable :: SPI(:,:,:), SPEI(:,:,:) ! (Nst,Ntime,Nlen) Real, Allocatable :: SPI_fcs(:,:,:), SPEI_fcs(:,:,:) Real, Allocatable :: Dr_Fcs(:,:,:) ! temporary array Character(len=20), Allocatable :: stname(:) ! Nst Real, Allocatable :: Lon(:),Lat(:),Reg(:) Real, Allocatable :: Ymax(:,:), Ymin(:,:) ! Nst,Nlen - Only used for Predictand Character (Len=5) :: yvarn !---------------------------------------------------------------------------------------------- Real, Allocatable :: X_inp(:,:), Y_inp(:), Yest(:,:), Yhq(:,:) !-------------------------------------------------------------------------------------- ! Tinh cac EOF cua cac ClimIDX: Real*8, Allocatable :: xdat(:,:), ydat(:,:), xtb(:), sx(:) Real*8, Allocatable :: EOFs(:,:), PCs(:,:) Real*8, Allocatable :: Eigenvalues(:), ECumulate(:) Real, Allocatable :: Xtb_Clim(:), Sx_Clim(:), PCs_Clim(:,:), Eig_Clim(:) Real, Allocatable :: lamda(:,:), PC1(:,:) Real, Allocatable :: EOF_Clim(:,:) Integer :: Nt, Nx ! Kich thuoc thuc cua mang X Integer :: Mcount, Ncount, switch, Standard, PC_Stard, Nmode, K Real*8 :: Factor Character (Len=200) :: FnameTC, FnameIDX, Predictand !----------------------------------------------------------------------------- Integer :: ireg, ilead Integer :: Nlayer ! Num of layers including Input, output and hidden laysers Integer :: Nepoch ! Num of epochs Integer :: Mx, Mvars double precision :: learn_rate ! Learning_Rate Integer, Allocatable :: Layer_size(:) ! (Nlayer) Logical :: training_mode double precision, Allocatable :: tmp(:) Integer :: ilr, inode, itest, Ntest, max_node, icase ! max_node: Maximum values among Mxvar, Myvar and number of nodes in Layer_size !----------------------------------------------------------------------------- ! Ann_hai: Integer :: Nsample, Nepoc !----------------------------------------------------------------------------- Character (Len=255) :: Dir_ANN_Save, save_file, Inp_Dir, Out_Dir, Mask_File, Fname Character (Len=255) :: DirInEC,DirInObs,Dir_Cli,Dir_Corrected Namelist /ann_info/ Nlayer,learn_rate,Nepoch,training_mode,yvarn Namelist /layer_info/ Layer_size ! Layer_size(1) will be Mxvar, Layer_size(Nlayer) will be Myvar !-------------------------------------------------------------------------------------- Integer :: imon, iyr, iidx, Nst Real, Allocatable :: Corr(:,:),Prob(:,:) ! (Nidx,Nidx) - Correlation between Clim_IDXa Real, Allocatable :: Cor_spi(:,:,:),Prb_spi(:,:,:) ! (Nst,Nlen,Nidx) - Correlation between Clim index and Dr index Real, Allocatable :: Cor_spei(:,:,:),Prb_spei(:,:,:) ! (Nst,Nlen,Nidx) Real, Allocatable :: Cor_spi_eof(:,:,:),Prb_spi_eof(:,:,:) ! (Nst,Nlen,Nmode) Real, Allocatable :: Cor_spei_eof(:,:,:),Prb_spei_eof(:,:,:) ! !----------------------------------------------------------------------------- Open (1, file="Namelist_Dr_fcs_Ann_new.txt",Status="old") Read (1, ann_info) ! Nlayer,learn_rate,Nepoch Allocate (Layer_size(Nlayer)) Read (1, layer_info) ! Close (1) !-------------------------------------------------------------------------------------- FnameTC = "/home/tanpv/Progs/4BSc/Res/spi_spei_stn_1980-2020.nc" FnameIDX = "/work/projects/G09/SEACLID/ClimIDX_2022/Clim_IDX_2022.nc" Dir_ANN_Save = "/home/tanpv/Progs/4BSc/ANN_Save/" Dir_Corrected = "/home/tanpv/Progs/4BSc/Res/" Out_Dir = "/home/tanpv/Progs/4BSc/4Grads/" Call Read_ClimIDX (FnameIDX) Call Read_Dr_IDX (FnameTC) !-------------------------------------------------------------------------------------- ! PCA of Climate indices: ! Tam thoi chua tinh ! Call PCA_of_ClimIDX !print*,"PC Clim:" !print "(A15,20(A11,I1.1))"," ",("eof_",i,i=1,Nmode) !Do j=1,Nidx ! print "(A15,20F12.3)",trim(Clim_nm(j)),(PCs_Clim(j,k),k=1,Nmode) !Enddo !-------------------------------------------------------------------------------------- Nspl = Ntime Ncase = Nst*Nlen*Nlag ! So lan tinh = So tram x So timescale x So lan du bao (han du bao) Mx=Nidx ! So bien doc lap = So climIDXs (hoac so (eofs) Mvars=Mx Mxvar=Mx; Myvar=1 !----------------------------------------------------------------------------- inode = Maxval(Layer_size(1:Nlayer)); max_node = Max(inode,Mxvar,Myvar) Allocate (Ymax(Nst,Nlen), Ymin(Nst,Nlen)) if (training_mode) then Allocate (Mnet_layers(Nlayer)) Allocate (Mnet_weight(Ncase,Nlayer-1,max_node,max_node)) Allocate (Mnet_biases(Ncase,Nlayer-1,max_node)) else Fname=trim(Dir_ANN_Save)//"Sav_"//trim(yvarn)//"_Ann_new.txt" Call Read_Ann_info (Fname) endif !----------------------------------------------------------------------------- Allocate (Dr_Fcs(Nst,Ntime,Nlen)) ! Forecast results icase = 0 Do ist=1,Nst Do ilen=1,Nlen Do ilag=Lag(1),Lag(Nlag) icase = icase+1 print*,icase,"/",Ncase Call Make_data (SPI,ist,ilen,ilag) ! Yobs is replated by SPI or SPEI Call Calc_ANN_New (training_mode,ist,ilen,ilag,icase) Enddo Enddo Enddo if (training_mode) then Fname=trim(Dir_ANN_Save)//"Sav_"//trim(yvarn)//"_Ann_new.txt" Call Write_Ann_info (Fname) endif if (training_mode) then ! Output filename Fname=trim(Dir_Corrected)//trim(yvarn)//"_Ann_new_training.nc" else Fname=trim(Dir_Corrected)//trim(yvarn)//"_Ann_new_testing.nc" endif !Fname="/home/tanpv/Progs/4BSc/Res/Fcs_spi_spei_stn_1980-2020.nc" Call Write_Dr_fcs_to_NC (Fname) !-------------------------------------------------------------------------------------- CONTAINS !---------------------------------------------------------------------------------------------- ! Real, Allocatable :: ClimIDX(:,:) ! (Nidx,Ntime ! Real, Allocatable :: SPI(:,:,:), SPEI(:,:,:) ! (Nst,Ntime,Nlen) Subroutine Make_data (Yobs,ist,ilen,ilag) Integer :: ist,ilen,ilag,it,i,j,idelta Real, intent (in) :: Yobs(:,:,:) !Nst,Ntime,Nlen) ! ilag: So thang tre cua ClimIDX so voi thoi diem Dr(i) duoc tinh (duoc du bao) --> ilag la LT ! Neu thoi diem (thang) duoc du bao la i, thoi diem lam du bao la i-ilead ! O day tam lay So thang bi tre SL la idelta=2 thang ! Clim: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ... ! Dr : 1 2 3 4 5 6 ... ! + + x ----LT----o ---> voi idelta=2 thi Leadtime = 4 = ilag ! + + x -------LT----o ---> voi idelta=2 thi Leadtime = 5 = ilag idelta = 2 k = idelta+ilag Nsample1=Nspl ! Dung luong mau if (Allocated(X_inp)) Deallocate (X_inp) if (Allocated(Y_inp)) Deallocate (Y_inp) Allocate (X_inp(Nsample1,Mx),Y_inp(Nsample1)) Do i=1,Nspl Y_inp(i) = Yobs(ist,i,ilen) ! SPEI(ist,i,ilen) Do j=1,Nidx X_inp(i,j) = ClimIDX(j,Nmon-k+i) Enddo !print"(3I5,20F7.3)",ilag,i,Nmon-k+i,Y_inp(i),(X_inp(i,j),j=1,6) Enddo !pause Nsample1=Nspl End Subroutine Make_data !-------------------------------------------------------------------------------- Subroutine Calc_ANN_New (training_mode,ist,ilen,ilag,icase) Integer :: icase,ist,ilen,ilag,ispl,N_spl, i,j,k,icount Real :: XX(Ntime,Mxvar),YY(Ntime),xtmp(Mxvar), ytmp Double precision :: x1d (Mxvar), y1d (Myvar) Logical :: training_mode,OK Integer, Allocatable :: Check(:) if (Allocated (Check)) Deallocate (Check) Allocate (Check(Ntime)) Do i=1,Ntime YY(i)=Y_inp(i) Do j=1,Mxvar XX(i,j)=X_inp(i,j) Enddo Enddo ! icount=0; Check(:)=0 Do i=1,Ntime ytmp=YY(i); xtmp(1:Mxvar)=XX(i,1:Mxvar) OK=Valid_Ar1(xtmp,Mxvar,Miss).AND.(ytmp/=Miss) if (OK) then icount=icount+1; Check(i)=1 endif Enddo N_spl = icount if (Allocated (Out_Y)) Deallocate (Out_Y) if (Allocated (Inp_X)) Deallocate (Inp_X) Allocate (Inp_X(N_spl,Mxvar),Out_Y(N_spl,Myvar)) ! To be used by subroutine neural_train_network_file_TP ispl=0 Do i=1,Ntime if (Check(i)==1) then ispl=ispl+1 Out_Y(ispl,1) = YY(i) Inp_X(ispl,1:Mxvar) = XX(i,1:Mxvar) endif Enddo !-------------------------------------------------------- ! Normalize: (for subroutine neural_train_network_file_TP) Ymax(ist,ilen) = Maxval(Out_Y(1:N_spl,1)) Ymin(ist,ilen) = Minval(Out_Y(1:N_spl,1)) Do i=1,N_spl Out_Y(i,1) = (Out_Y(i,1)-Ymin(ist,ilen))/(Ymax(ist,ilen)-Ymin(ist,ilen)) Enddo !-------------------------------------------------------- Layer_size(1) = Mxvar; Layer_size(Nlayer) = Myvar ! Update Num of nodes for input (X) and output (Y) layers call neural_init_network(Net, Layer_size) if (training_mode) then call neural_train_network_file_TP(Net, Nepoch, learn_rate) call neural_save_network_TP (icase, net) else call neural_load_network_TP(icase,Nlayer,net) endif if (Allocated (Yest)) Deallocate (Yest) Allocate (Yest(Ntime,Myvar)); Yest = Miss ispl=0 Do i=1,Ntime if (Check(i)==1) then ispl=ispl+1 x1d(1:Mxvar) = Inp_X(ispl,1:Mxvar) y1d(:) = neural_use_network(Net,x1d) Do j=1,Myvar Yest(i,j) = y1d(j)*(Ymax(ist,ilen)-Ymin(ist,ilen)) + Ymin(ist,ilen) Enddo !print"(A,6F8.4)","Yest,Y_inp: ",Yest(i,:), YY(i) endif Enddo !EC_bc(ilon,ilat,1:Nspl,ilead) = Yest(1:Ntime,1) Dr_Fcs(ist,1:Ntime,ilen) = Yest(1:Ntime,1) End Subroutine Calc_ANN_New !-------------------------------------------------------------------------------- Subroutine Write_Ann_info (Fname) Character (Len=*) :: Fname Open (1,file=Fname,form="unformatted", status="unknown") write (1) Nlayer, max_node, Ncase write (1) Mnet_layers write (1) Mnet_weight write (1) Mnet_biases End Subroutine Write_Ann_info !-------------------------------------------------------------------------------- Subroutine Read_Ann_info (Fname) Character (Len=*) :: Fname Open (1,file=Fname,form="unformatted", status="old") read (1) Nlayer, max_node, Ncase Allocate (Mnet_layers(Nlayer)) Allocate (Mnet_weight(Ncase,Nlayer-1,max_node,max_node)) Allocate (Mnet_biases(Ncase,Nlayer-1,max_node)) read (1) Mnet_layers read (1) Mnet_weight read (1) Mnet_biases End Subroutine Read_Ann_info !-------------------------------------------------------------------------------------- Subroutine PCA_of_ClimIDX Integer :: it, imode,iidx, i Real*8 :: Factor switch=0; Standard=1; PC_Stard=0; Nmode=6 Nt=Ntime1; Nx=Nidx 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)) Do it=1,Nt !=Ntime Do i=1,Nx ! =Nidx*Nlag xdat(it,i)=ClimIDX(i,it) Enddo Enddo Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, & Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode) print"(A,6(I2,F6.2),F6.2)","Eigen - mode ",(imod, Eigenvalues(imod)/sum(Eigenvalues)*100,imod=1,Nmode),sum(Eigenvalues(1:Nmode))/sum(Eigenvalues)*100. Do imod=1,Nmode Do it=1,Nt EOF_Clim(it,imod)=EOFs(it,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 End Subroutine PCA_of_ClimIDX !-------------------------------------------------------------------------------------- Subroutine write_Correlation_Dr_eof_of_Clim_to_GrADS (Cor,Prb,Fname) Character (Len=*) :: Fname Real, Dimension (Nst,Nlen,Nmode) :: Cor,Prb Real :: X(Nst,Nidx,Nlen*2) Integer :: ist,idx,Nlen2 ! Nlen2=Nlen*2 ! Do ist=1,Nst ! Do idx=1,Nmode ! X(ist,idx,1) = Cor(ist,1,idx) ! X(ist,idx,2) = Cor(ist,2,idx) ! X(ist,idx,3) = Prb(ist,1,idx) ! X(ist,idx,4) = Prb(ist,2,idx) ! Enddo ! Enddo ! Call Write_GrADS_Stn_Multivar (X,Nst,Nmode,Nlen2,Lon,Lat,stname,"spi",Fname) End Subroutine write_Correlation_Dr_eof_of_Clim_to_GrADS !-------------------------------------------------------------------------------------- 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,Ntime1)) 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=Yr1-1,Yr2 ! Voi ClimIDX thi Timeseries dai hon 12 thang truoc timeseries cua Drought 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 Deallocate (Year,Month,tmp) Enddo End Subroutine Read_ClimIDX !-------------------------------------------------------------------------------------- Subroutine Read_Dr_IDX (Fname) Character (Len=*) :: Fname Real, allocatable :: stnid(:),time(:),xlen(:),X3d(:,:,:) Integer :: Nt,Nle,it1 !Allocate (SPI(Nst,Ntime,Nlen),SPEI(Nst,Ntime,Nlen)) !SPI = Miss; SPEI = Miss it1 = Nmon+1 ! Tu thang thu 13 (bo nam dau tien) Call rd_ncdf_3d (Fname,stnid,time,xlen,Nst,Nt,Nle,X3d,"spi") Allocate (SPI(Nst,Ntime,Nlen),SPEI(Nst,Ntime,Nlen)) SPI = Miss; SPEI = Miss Where (IsNaN(X3d)) X3d=Miss SPI(:,1:Ntime,:) = X3d(:,it1:Nt,:) Deallocate (stnid,time,xlen,X3d) Call rd_ncdf_3d (Fname,stnid,time,xlen,Nst,Nt,Nle,X3d,"spei") Where (IsNaN(X3d)) X3d=Miss SPEI(:,1:Ntime,:) = X3d(:,it1:Nt,:) Deallocate (stnid,time,xlen,X3d) Call rd_ncdf_1d (Fname,stnid,Nst,Lon,"lon") Deallocate (stnid) Call rd_ncdf_1d (Fname,stnid,Nst,Lat,"lat") Deallocate (stnid) Call rd_ncdf_1d (Fname,stnid,Nst,Reg,"reg") Deallocate (stnid) Call rd_ncdf_1d_Char (Fname,stnid,Nst,stname,"stname") Deallocate (stnid) End Subroutine Read_Dr_IDX !------------------------------------------------------------------------------ Subroutine Write_Dr_fcs_to_NC (Fname) Character (Len=*) :: Fname Real, Allocatable :: Stn(:), time(:), LenDr(:), X2D(:,:) Integer :: it, Nm, ist, imon, iyr, Nstn Character (Len=255) :: var_nm,t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm Allocate (Stn(Nst),time(Ntime),LenDr(NLen)) Do ist=1,Nst; Stn(ist)=ist; Enddo Do it=1,Ntime; time(it)=it; Enddo Do it=1,NLen; LenDr(it)=Len_SPI(it); Enddo ! SPI t_unit="months since 1980-01-01 "; lon_nm="stn"; lat_nm="mon"; t_nm="len"; var_unit="-" var_ln="SPI for period 1980-2020 with drought timescale=1,3,6,12 months" Call wr_ncdf_3d (Fname,Stn,time,LenDr,Nst,Ntime,Nlen,Dr_Fcs,"spi",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! SPEI !var_ln="SPEI for period 1980-2020 with drought timescale=1,3,6,12 months, PET-thornthwaite" !Call wr_ncdf_3d_apd (Fname,SPEI_th,"spei",Nst,Ntime,Nlen, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! Lon, Lat, Reg and Station name var_ln="Longitude in East degree"; var_unit="deg-E" Call wr_ncdf_1d_apd (Fname,Lon,"lon", Nst, var_ln, var_unit, lon_nm, Miss) var_ln="Latitude in North degree"; var_unit="deg-N" Call wr_ncdf_1d_apd (Fname,Lat,"lat", Nst, var_ln, var_unit, lon_nm, Miss) var_ln="Climatic sub-regions"; var_unit="-" Call wr_ncdf_1d_apd (Fname,Reg,"reg", Nst, var_ln, var_unit, lon_nm, Miss) var_ln="Name of Stations"; var_unit="st" Call wr_ncdf_1d_apd_char (Fname,Stname,"stname", Nst, var_ln, var_unit, lon_nm, Miss) End Subroutine Write_Dr_fcs_to_NC !-------------------------------------------------------------------------------------- Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Function_Valid_Ar1.f90" Include "/home/tanpv/Mylibs/Calc_Correlation_Coefficient.f90" Include "/home/tanpv/Mylibs/Sub_PCA_2019.f90" Include "/home/tanpv/Mylibs/Write_GrADS_Stn_Multivar.f90" !-------------------------------------------------------------------------------------- END !-------------------------------------------------------------------------------------- Include "/home/tanpv/Mylibs/Student_Prob.f90" !--------------------------------------------------------------------------------------