! ! CHUONG TRINH thu nghiem hieu chinh KQ mo hinh Hybrid Statistical-Dynamical method ! Cac Mo hinh thong ke duoc su dung: ! 1) ANN: ! 2) HQTT: Hoi quy tuyen tinh nhieu bien (Ket qua giong voi HQTB) ! 3) LAD: Least Absolute Deviation ! 4) LMV: Minimax ! Cac nhan to du bao: Cac bien du bao cua mo hinh PROGRAM BC_ANN Use Netcdf Integer, Parameter :: Yr1=2011, Yr2=2021, Nmon=12, Nreg=7, NLead=6, Nvar=4, Nfile=9 Integer, Parameter :: Nyrs=(Yr2-Yr1+1), Ntime=Nyrs*Nmon, Nmem=6, Nmem1=7 Real, Parameter :: Miss=-99. Character (Len=3) :: vname(Nvar)=(/"obs","ann","mod","cli"/) Character (Len=255) :: Fname Real, Allocatable :: lon(:),lat(:),Regs(:,:), X3d(:,:,:) Real, Allocatable :: xlon(:),xlat(:),xlev(:),xtime(:) Integer :: Nlo,Nla,Nz,Nt Integer :: Nlon, Nlat, Ngrd, ilon, ilat, it Integer :: Mvars Real, Allocatable :: Dat_All(:,:,:,:,:,:) ! (Nlon,Nlat,Ntime,Nvar,Nmem1,NLead) Real, Allocatable :: X_inp(:,:), Y_inp(:), Yest(:,:), Yhq(:,:) Character (Len=5) :: yvarn, Model_nm !-------------------------------------------------------------------------------- ! Verification: Real, Allocatable, Dimension(:,:,:,:,:) :: Pr_ME, Pr_MAE, Pr_RME, Pr_RMAE, Pr_COR, Pr_Pval, Pr_AV ! (Nlon,Nlat,Nmon,NLead,Nvar) Integer :: Nyear, imon ! Tong so nam co SL va chi so thang !-------------------------------------------------------------------------------- ! ANN: character (len=250) :: input_file, output_file, save_file character (len=250) :: Dir, Dir_ANN_Out, Dir_Cli, Dir_ANN_Save, Dir_Veri, Dir_4Dr logical :: training_mode ! Use for calculating scaling factors Integer :: Nsample, Nepoc, Nsample1 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,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_rate0, learning_rate, momentum, error, tmp, mintmp logical :: debug, Run_Cli !-------------------------------------------------------------------------------- Integer :: Op_Nhidden, Op_Nepoc ! Optimal values Real :: Op_learning_rate Real, Allocatable :: X4d(:,:,:,:), X5d(:,:,:,:,:), F1d(:),O1d(:), Score_Cor(:,:,:,:), Score_AV(:,:,:,:) !------------------- For ANN ------------------------------------------------------------- Namelist /ann_info/ Nhidden,learning_rate,momentum,training_mode,Nepoc,debug Namelist /Model_info/ Model_nm, yvarn ! Nhidden = 5 !4 ! Number of hidden layers ! learning_rate = 0.01 ! momentum = 0.9 ! 10.0 ! training_mode = .true. ! For training ! Nepoc = 1000 ! Iterative training times ! debug = .false. !-------------------------------------------------------------------------------- Open (1, file="Namelist_ANN_Training.txt",Status="old") Read (1, ann_info) Read (1, Model_info) learning_rate0 = learning_rate Fname="/home/tanpv/Data/AWC/Mask_VN_0.25_Corrected.nc" Call rd_ncdf_2d (Fname,lon,lat,Nlon,Nlat,Regs,"reg") ! 7 sub-regs mask !Model_nm="clWRF"; yvarn="pr" !Dir_ANN_Out ="/home/tanpv/Data/562_VNU/0.25_VN/ANN_Out/"; !Dir_ANN_Save="/home/tanpv/Data/562_VNU/0.25_VN/ANN_Save/" !Dir_Cli ="/work/projects/G09/562_VNU/0.25_VN/Cli_Out/" !"/home/tanpv/Data/562_VNU/0.25_VN/Cli_Out/" Dir_Veri ="/home/tanpv/Data/562_VNU/0.25_VN/Verify/" Dir_4Dr ="/home/tanpv/Data/562_VNU/0.25_VN/4Drou/" !-------------------------------------------------------------------------------- Call Read_all_data (Dir_4Dr,Model_nm,yvarn) Nvar1=Nvar-1 ! ann, model, cli !-------------------------------------------------------------------------------- Allocate (Pr_RME(Nlon,Nlat,NLead,Nmon,Nvar1),Pr_RMAE(Nlon,Nlat,NLead,Nmon,Nvar1),Pr_COR(Nlon,Nlat,NLead,Nmon,Nvar1),Pr_Pval(Nlon,Nlat,NLead,Nmon,Nvar1)) Allocate (Pr_ME (Nlon,Nlat,NLead,Nmon,Nvar1),Pr_MAE (Nlon,Nlat,NLead,Nmon,Nvar1),Pr_AV (Nlon,Nlat,NLead,Nmon,Nvar1)) Pr_ME=Miss; Pr_MAE=Miss; Pr_RME=Miss; Pr_RMAE=Miss; Pr_COR=Miss; Pr_Pval=Miss; Pr_AV=Miss Allocate (Score_Cor(Nreg,Nmon,NLead,Nvar1), Score_AV(Nreg,Nmon,NLead,Nvar1)) ! Number of gridpoints with Cor>0, AV>0 Nyear=Nyrs Do imem=1,Nmem1 Call Verification(imem) write(Fname,"(A,I2.2,A)") trim(Dir_Veri)//"Verification_scores_testing_"//trim(yvarn)//"_"//trim(Model_nm)//"_mem_",imem,".nc" Call Write_Verification (Fname,yvarn) write(Fname,"(A,I2.2,A)") trim(Dir_Veri)//"Score_Cor_AV_testing_"//trim(yvarn)//"_"//trim(Model_nm)//"_mem_",imem,".txt" Open (3,file=Fname,status="Unknown") Do ivar=1,Nvar1 write (3,"(A)") "Percentage of grid points with positive correlation and added value of "//trim(vname(ivar+1)) Do ilead=1,NLead write (3,"(A,I1.1,A)") "Lead time = ",ilead," month(s)" write (3,"(6x,A6,12I6)") " Reg. ",(imon, imon=1,Nmon) Do ireg=1,Nreg write (3,"(A6,I6,12F6.1)") "Cor - ",ireg,(Score_Cor(ireg,imon,ilead,ivar),imon=1,Nmon) write (3,"(2A6,12F6.1)") "AV - "," ",(Score_AV (ireg,imon,ilead,ivar),imon=1,Nmon) Enddo Enddo Enddo Close (3) Enddo !-------------------------------------------------------------------------------- CONTAINS !-------------------------------------------------------------------------------- Subroutine Verification (imem) Integer :: ilon,ilat,imon,imon1,ilead,ireg, ic, ic0, ic1, ic2, it Integer, Parameter :: N=Nyrs ! Dat_All(Nlon,Nlat,Ntime,Nvar,Nmem1,NLead) Real :: obs(N),ann(N),model(N),cli(N),Corr,Prob Do ilead=1,NLead Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))>0) then Do imon=1,Nmon ic=0 Do it=1,Ntime imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon if (imon1==imon) then ic=ic+1 obs(ic) =Dat_All(ilon,ilat,it,1,imem,ilead); ann(ic)=Dat_All(ilon,ilat,it,2,imem,ilead) model(ic)=Dat_All(ilon,ilat,it,3,imem,ilead); cli(ic)=Dat_All(ilon,ilat,it,4,imem,ilead) endif Enddo if (ic>N) then; print*,"Something's wrong in Verification"; stop; endif Pr_ME (ilon,ilat,ilead,imon,1) = ME (ann, obs,ic,Miss) Pr_ME (ilon,ilat,ilead,imon,2) = ME (model,obs,ic,Miss) Pr_ME (ilon,ilat,ilead,imon,3) = ME (cli, obs,ic,Miss) Pr_MAE (ilon,ilat,ilead,imon,1) = MAE (ann, obs,ic,Miss) Pr_MAE (ilon,ilat,ilead,imon,2) = MAE (model,obs,ic,Miss) Pr_MAE (ilon,ilat,ilead,imon,3) = MAE (cli, obs,ic,Miss) Pr_RME (ilon,ilat,ilead,imon,1) = RME_1 (ann, obs,ic,Miss) Pr_RME (ilon,ilat,ilead,imon,2) = RME_1 (model,obs,ic,Miss) Pr_RME (ilon,ilat,ilead,imon,3) = RME_1 (cli, obs,ic,Miss) Pr_RMAE(ilon,ilat,ilead,imon,1) = RMAE_1(ann, obs,ic,Miss) Pr_RMAE(ilon,ilat,ilead,imon,2) = RMAE_1(model,obs,ic,Miss) Pr_RMAE(ilon,ilat,ilead,imon,3) = RMAE_1(cli, obs,ic,Miss) Call Correlation (ann, obs,ic,Corr,Prob,Miss); Pr_COR(ilon,ilat,ilead,imon,1)=Corr; Pr_Pval(ilon,ilat,ilead,imon,1)=Prob Call Correlation (model,obs,ic,Corr,Prob,Miss); Pr_COR(ilon,ilat,ilead,imon,2)=Corr; Pr_Pval(ilon,ilat,ilead,imon,2)=Prob Call Correlation (cli, obs,ic,Corr,Prob,Miss); Pr_COR(ilon,ilat,ilead,imon,3)=Corr; Pr_Pval(ilon,ilat,ilead,imon,3)=Prob Pr_AV (ilon,ilat,ilead,imon,1) = AV (ann,model,obs,ic,Miss) ! Added value of ann vs model Pr_AV (ilon,ilat,ilead,imon,2) = AV (cli,model,obs,ic,Miss) ! Added value of cli vs model Pr_AV (ilon,ilat,ilead,imon,3) = Miss ! Not use Enddo endif Enddo Enddo Enddo ! Do ivar=1,Nvar1 Do ilead=1,NLead Do imon=1,Nmon Do ireg=1,Nreg ic0=0; ic1=0; ic2=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))==ireg) then ic0=ic0+1 if (Pr_COR (ilon,ilat,ilead,imon,ivar)>0.2) ic1=ic1+1 if (Pr_AV (ilon,ilat,ilead,imon,ivar)>0.1) ic2=ic2+1 endif Enddo Enddo Score_Cor(ireg,imon,ilead,ivar)=real(ic1)/real(ic0)*100.0 if (ivar