! 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) ! double precision, Allocatable :: tmp(:) Integer :: ilr, inode, itest, Ntest, max_node ! max_node: Maximum values among Mxvar, Myvar and number of nodes in Layer_size !----------------------------------------------------------------------------- ! Ann_hai: Integer :: Nsample, Nepoc character (len=250) :: input_file, Out_file, save_file logical :: training_mode Real,allocatable, Dimension(:,:) :: Inputs, Outputs, Inputs_X, Outputs_Y Real,allocatable, Dimension(:) :: min_in, max_in, min_out, max_out Real :: mean_error, mean_error_check integer :: Nin, Nhidden, Nout, iin, ihid, iout integer :: i,j 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_rate0, learning_rate,momentum, error, tmp, mintmp logical :: debug !------------------------------------------------------------------------------- ! controller_save_info: ! Nin,Nhidden,Nout - Nin=MXvar; Nout=MYvar Integer :: Ncase,icase Real,allocatable, Dimension(:,:,:) :: weight_ip ! (Ncase,Nin,Nhidden) - Ncase=So lan tinh, Real,allocatable, Dimension(:,:,:) :: weight_op ! (Ncase,Nhidden,Nout) Real,allocatable, Dimension(:,:) :: min_ip, max_ip, min_op, max_op ! (Ncase) ! For saving in ANN subroutine: Real,allocatable, Dimension(:,:) :: wei_i, wei_o Real,allocatable, Dimension(:) :: min_i, max_i, min_o, max_o !----------------------------------------------------------------------------- Character (Len=255) :: Dir_ANN_Save, Inp_Dir, Out_Dir, Mask_File, Fname Character (Len=255) :: DirInEC,DirInObs,Dir_Cli,Dir_Corrected Namelist /ann_info/ Nhidden,learning_rate,momentum,training_mode,Nepoc,debug Namelist /Model_info/ yvarn !----------------------------------------------------------------------------- Open (1, file="Namelist_Dr_fcs_Ann_Hai.txt",Status="old") Read (1, ann_info) Read (1, Model_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) !----------------------------------------------------------------------------- if (training_mode) then MXvar = Nidx; MYvar = 1 Mvars = MXvar Mx=Nidx Allocate (weight_ip(1:Ncase,0:MXvar,1:Nhidden),weight_op(Ncase,0:Nhidden,MYvar)) Allocate (min_ip(Ncase,MXvar),max_ip(Ncase,MXvar),min_op(Ncase,MYvar),max_op(Ncase,MYvar)) else Fname=trim(Dir_ANN_Save)//"Sav_"//trim(yvarn)//".txt" Call Read_Ann_info (Fname) Mvars = MXvar endif Allocate (wei_i(0:MXvar,1:Nhidden), wei_o(0:Nhidden,MYvar)) ! To be used in ANN Sub Allocate (min_i(MXvar),max_i(MXvar),min_o(MYvar),max_o(MYvar)) ! To be used in ANN Sub !----------------------------------------------------------------------------- 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) Call ANN_2025(Inputs_X,Outputs_Y,Nsample1,Mvars,1,Miss) !ANN_2023 (Xinp,Youtp,Nspl,MXvar,MYvar,Miss) if (training_mode) then Call Save_Ann_info (icase) else Call Load_Ann_info (icase) endif Dr_Fcs(ist,1:Ntime,ilen) = Yest(1:Ntime,1) ! EC_bc(ilon,ilat,1:Nspl,ilead) = Yest(1:Nspl,1) ! Ann outputs Enddo Enddo Enddo if (training_mode) then Fname=trim(Dir_ANN_Save)//"Sav_"//trim(yvarn)//"_Ann_Hai.txt" Call Write_Ann_info (Fname) endif if (training_mode) then ! Output filename Fname=trim(Dir_Corrected)//trim(yvarn)//"_Ann_Hai_training.nc" else Fname=trim(Dir_Corrected)//trim(yvarn)//"_Ann_Hai_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 ! Neu thoi diem (thang) duoc du bao la i, thoi diem lam du bao la i-ilead thi ilag=i-ilead-So thang bi tre SL ! O day tam lay So thang bi tre SL la idelta=2 thang 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 if (Allocated(Inputs_X)) Deallocate (Inputs_X) if (Allocated(Outputs_Y)) Deallocate (Outputs_Y) Allocate (Inputs_X(Mx,Nsample1), Outputs_Y(1,Nsample1)) i=0 Do it=1,Nspl i=i+1 Outputs_Y(1,i)=Y_inp(it) Do j=1,Mx Inputs_X(j,i)=X_inp(it,j) Enddo Enddo End Subroutine Make_data !-------------------------------------------------------------------------------- Subroutine Save_Ann_info (icase) Integer :: icase ! Nin=MXvar; Nout=MYvar wei_i, wei_o weight_ip(icase,0:MXvar,1:Nhidden) = wei_i(0:MXvar,1:Nhidden) weight_op(icase,0:Nhidden,1:MYvar) = wei_o(0:Nhidden,1:MYvar) min_ip(icase,1:MXvar) = min_i(1:MXvar); max_ip(icase,1:MXvar) = max_i(1:MXvar) min_op(icase,1:MYvar) = min_o(1:MYvar); max_op(icase,1:MYvar) = max_o(1:MYvar) End Subroutine Save_Ann_info !---------------------------------------------------------------------------------------------- Subroutine Load_Ann_info (icase) Integer :: icase ! Nin=MXvar; Nout=MYvar wei_i, wei_o wei_i(0:MXvar,1:Nhidden) = weight_ip(icase,0:MXvar,1:Nhidden) wei_o(0:Nhidden,1:MYvar) = weight_op(icase,0:Nhidden,1:MYvar) min_i(1:MXvar) = min_ip(icase,1:MXvar); max_i(1:MXvar) = max_ip(icase,1:MXvar) min_o(1:MYvar) = min_op(icase,1:MYvar); max_o(1:MYvar) = max_op(icase,1:MYvar) End Subroutine Load_Ann_info !-------------------------------------------------------------------------------- Subroutine Read_Ann_info (Fname) Character (Len=*) :: Fname Open (1,file=trim(Fname),form="unformatted",status="old") read(1) MXvar,Nhidden,MYvar,Ncase Allocate (weight_ip(1:Ncase,0:MXvar,1:Nhidden),weight_op(Ncase,0:Nhidden,MYvar)) Allocate (min_ip(Ncase,MXvar),max_ip(Ncase,MXvar),min_op(Ncase,MYvar),max_op(Ncase,MYvar)) read(1) weight_ip read(1) weight_op read(1) min_ip read(1) max_ip read(1) min_op read(1) max_op Close (1) End Subroutine Read_Ann_info !---------------------------------------------------------------------------------------------- Subroutine Write_Ann_info (Fname) Character (Len=*) :: Fname Open (1,file=trim(Fname),form="unformatted",status="unknown") write(1) MXvar,Nhidden,MYvar,Ncase write(1) weight_ip write(1) weight_op write(1) min_ip write(1) max_ip write(1) min_op write(1) max_op Close (1) End Subroutine Write_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/Sub_neural_2025.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" !--------------------------------------------------------------------------------------