! ! 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=1982, Yr2=2011, Nmon=12, Nreg=7, NLead=5, Nvar=17, Nfile=9 Integer, Parameter :: Nyrs=(Yr2-Yr1+1) Real, Parameter :: Miss=-99. Character (Len=4) :: vname(Nvar)=(/"pr_o","t2_o","tx_o","tm_o","pr","t2","tx","tm", & "h500","h850","olr","q2","rh2","slp","u200","u850","w700"/) Integer :: Fcs_mon, Nyr1, NMax, Ncase, Ntc_var Character (Len=255) :: Fname Real, Allocatable :: lon(:),lat(:),Regs(:,:), X3d(:,:,:) ! xlon,xlat,xlev,xtime,Nlo,Nla,Nz,Nt Real, Allocatable :: xlon(:),xlat(:),xlev(:),xtime(:) Integer :: Nlo,Nla,Nz,Nt Integer :: Nlon, Nlat, Ntime, Ngrd, ilon, ilat, it Integer :: Mvars Real, Allocatable :: Dat_All (:,:,:,:,:) ! (Nlon,Nlat,Ntime,Nvar,NLead) Real, Allocatable :: Dat_All_0(:,:,:,:,:) ! (Nlon,Nlat,Ntime,Nvar,NLead) Real, Allocatable :: Dat_mean(:,:,:,:,:) ! Dat_mean(Nlon,Nlat,Nmon,Nvar,NLead),Dat_std(Nlon,Nlat,Nmon,Nvar,NLead) Real, Allocatable :: Dat_std (:,:,:,:,:) Real, Allocatable :: Del_Cli(:,:,:,:) ! (Nlon,Nlat,NLead,Nmon) : Del = F - O Real, Allocatable :: Bias_Cli(:,:,:,:) ! (Nlon,Nlat,NLead,Nmon) - Model Cli. Bias = Model_Fcs - Obs ==> Corrected=Model_Fcs-Bias 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,NLead,Nmon) 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="/home/tanpv/Progs/562_VNU/Data/For_BC/" Call Read_all_data (Dir,Model_nm) 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/" !-------------------------------------------------------------------------------- ! Tinh Model Climatology: Run_Cli=.true. !Run_Cli=.false. if (Run_Cli) then Call Calc_Bias_Cli (yvarn) write (Fname,"(3(A,I2.2),A)") trim(Dir_Cli)//"Bias_Cli_"//trim(Model_nm)//"_"//trim(yvarn)//".nc" Call Write_Bias_Cli (Fname) if (Allocated(Bias_Cli)) Deallocate(Bias_Cli) Allocate (Bias_Cli(Nlon,Nlat,NLead+2,Nmon)); Bias_Cli=Del_Cli ! NLead+2: Mean Obs !stop else Fname=trim(Dir_Cli)//"Bias_Cli_"//trim(Model_nm)//"_"//trim(yvarn)//".nc" Call rd_ncdf_4d (Fname,xlon,xlat,xlev,xtime,Nlo,Nla,Nz,Nt,Bias_Cli,"bias") ! Model Cli. Bias = Model_Fcs - Obs ==> Corrected=Model_Fcs-Bias Deallocate (xlon,xlat,xlev,xtime) endif !-------------------------------------------------------------------------------- ! Data stadardization: print*,"Beginning Data standardization..." Call Data_standardization Fname=trim(Dir_Cli)//"Mean_and_Std_"//trim(Model_nm)//"_"//trim(yvarn)//".nc" Call Write_Mean_Std (Fname) !-------------------------------------------------------------------------------- ! Allocate (Pr_RME(Nlon,Nlat,NLead,Nmon),Pr_RMAE(Nlon,Nlat,NLead,Nmon),Pr_COR(Nlon,Nlat,NLead,Nmon),Pr_Pval(Nlon,Nlat,NLead,Nmon)) ! Allocate (Pr_ME(Nlon,Nlat,NLead,Nmon),Pr_MAE(Nlon,Nlat,NLead,Nmon),Pr_AV(Nlon,Nlat,NLead,Nmon)) ! 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), Score_AV(Nreg,Nmon,NLead)) ! Number of gridpoints with Cor>0, AV>0 Nyear=Nyrs Do ilead=1,NLead if (Allocated(X5d)) Deallocate (X5d) Allocate (X5d(Nlon,Nlat,Nmon,Nyear,4)) ! 5th-dimension = 4: 1) OBS, 2) ANN est, 3) Cli BC, 4) Model out X5d=Miss Do imon=1,Nmon if (Allocated(X4d)) Deallocate (X4d) Allocate (X4d(Nlon,Nlat,3,Nyrs)) !-1)) ! z=3: Yobs, Yest, Ymodel X4d=Miss !write (output_file,"(2(A,I2.2),A)") trim(Dir_ANN_Out)//"Out_Training_"//trim(yvarn)//"_"//trim(Model_nm)//"_mon_",imon,"_lead_",ilead,".nc" Do ireg=1,Nreg Call Make_Data (trim(yvarn),ireg,imon,ilead) write (save_file,"(3(A,I2.2),A)") trim(Dir_ANN_Save)//"Sav_"//trim(yvarn)//"_"//trim(Model_nm)//"_Reg_",ireg,"_mon_",imon,"_lead_",ilead,".txt" learning_rate = learning_rate0 ! 0.05 Call ANN_2023(Inputs_X,Outputs_Y,Nsample1,Mvars,1,Miss) !ANN_2023 (Xinp,Youtp,Nspl,MXvar,MYvar,Miss) Call Convert_to_lat_lon (trim(yvarn),ireg,imon,ilead) ! (ireg) Enddo Call Store_5th_Arr (ilead,imon) Enddo write (Fname,"(A,I2.2,A)") trim(Dir_4Dr)//"Dr_Training_"//trim(yvarn)//"_"//trim(Model_nm)//"_lead_",ilead,".nc" Call Write_Out_4_Drought (Fname,yvarn) Enddo !-------------------------------------------------------------------------------- CONTAINS !-------------------------------------------------------------------------------- Subroutine Write_Out_4_Drought (Fname,yvarn) Character (Len=*) :: Fname,yvarn Integer, parameter :: Nvar=4 Character (Len=250) :: var_nm, t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, var_longname Integer :: imon,iyr,imem,ivar,it, Nz Real, Allocatable :: Y4d(:,:,:,:), zlev(:), xtime(:) Character (Len=3), parameter :: vname(Nvar)=(/"obs","ann","mod","cli"/) Character (Len=20), parameter :: v_ln(Nvar)=(/"observation","ANN estimation","Model Forecast","Climatology bias correction"/) ! X6d(Nlon,Nlat,Nmon,Nyear,Nmem+1,4) Nz=Nvar; Allocate (zlev(Nz),xtime(Ntime),Y4d(Nlon,Nlat,Nz,Ntime)) Do i=1,Nz; zlev(i)=i; Enddo; Do i=1,Ntime; xtime(i)=i; Enddo t_unit="months since 1982-01-01 "; lon_nm="lon"; lat_nm="lat"; t_nm="month"; z_nm="member" var_ln="Outputs of "//trim(yvarn)//" for drought calculation. z=1:obs, 2:ann, 3:model, 4:cli" Do ivar=1,Nvar it=0 Do iyr=1,Nyear Do imon=1,Nmon it=it+1; Do ilon=1,Nlon Do ilat=1,Nlat Y4d(ilon,ilat,ivar,it)=X5d(ilon,ilat,imon,iyr,ivar) Enddo Enddo Enddo Enddo Enddo Call wr_ncdf_4d (Fname,Lon,Lat,zlev,xtime,Nlon,Nlat,Nz,Ntime,Y4d,trim(yvarn), & t_unit, var_ln, "mm/month", lon_nm, lat_nm, z_nm, t_nm, Miss) End Subroutine Write_Out_4_Drought !-------------------------------------------------------------------------------- Subroutine Write_Bias_Cli (Fname) Character (Len=*) :: Fname Character (Len=250) :: vname, t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm Real, Allocatable, Dimension (:) :: zlev, xtime INteger :: i, Nz, NLead2 NLead2=NLead+2 Allocate (zlev(NLead2),xtime(Nmon)) Do i=1,NLead2; zlev(i)=i; Enddo; Do i=1,Nmon; xtime(i)=i; Enddo t_unit="months since 1982-01-01 "; lon_nm="lon"; lat_nm="lat"; t_nm="mon"; z_nm="lead" var_ln="Differences between model outputs and observarions used as climatological biases. z=6 same as z=5, z=7: Mean Obs." Call wr_ncdf_4d (Fname,Lon,Lat,zlev,xtime,Nlon,Nlat,NLead2,Nmon,Del_Cli,"bias", & t_unit, var_ln, "-", lon_nm, lat_nm, z_nm, t_nm, Miss) print*,"Write_Bias_Cli... Done!" End Subroutine Write_Bias_Cli !-------------------------------------------------------------------------------- Subroutine Write_Mean_Std (Fname) Character (Len=*) :: Fname Character (Len=250) :: vname, t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm Character (Len=2), Parameter :: mon(Nmon)=(/"01","02","03","04","05","06","07","08","09","10","11","12"/) Real, Allocatable, Dimension (:) :: zvar, xtime Real, Allocatable, Dimension (:,:,:,:) :: A4d INteger :: imon,i, Nz, NLead2 NLead2=NLead+2 Allocate (zvar(Nvar),xtime(NLead2), A4d(Nlon,Nlat,Nvar,NLead2)) Do i=1,Nvar; zvar(i)=i; Enddo; Do i=1,NLead2; xtime(i)=i; Enddo t_unit="months since 1982-01-01 "; lon_nm="lon"; lat_nm="lat"; t_nm="lead"; z_nm="var" Do imon=1,Nmon var_ln="Mean values over entire baseline timeseris (1982-2011). Toal LT is 5, LT=6,7 are signed by LT=5. Month "//mon(imon) A4d(:,:,:,:)=Dat_mean(:,:,imon,:,:); vname="mean"//mon(imon) if (imon==1) then Call wr_ncdf_4d (Fname,Lon,Lat,zvar,xtime,Nlon,Nlat,Nvar,NLead2,A4d,trim(vname), & t_unit, var_ln, "-", lon_nm, lat_nm, z_nm, t_nm, Miss) else Call wr_ncdf_4d_apd (Fname,A4d,trim(vname),Nlon,Nlat,Nvar,NLead2, & var_ln, "-", lon_nm, lat_nm, z_nm, t_nm, Miss) endif Enddo Do imon=1,Nmon var_ln="Stadard deviation values over entire baseline timeseris (1982-2011). Toal LT is 5, LT=6,7 are signed by LT=5. Month "//mon(imon) A4d(:,:,:,:)=Dat_std(:,:,imon,:,:); vname="std"//mon(imon) Call wr_ncdf_4d_apd (Fname,A4d,trim(vname),Nlon,Nlat,Nvar,NLead2, & var_ln, "-", lon_nm, lat_nm, z_nm, t_nm, Miss) Enddo Deallocate (zvar, xtime, A4d) End Subroutine Write_Mean_Std !-------------------------------------------------------------------------------- Subroutine Convert_to_lat_lon (yvarn,ireg,imon,ilead) !(ireg) Character (Len=*) :: yvarn Integer :: ireg,imon,ilead,ilon,ilat,it,iyr,iobs,imod Real :: Omax, Omin, Fmax, Fmin if (yvarn=="pr") then iobs=1; imod=5 elseif (yvarn=="t2") then iobs=2; imod=6; elseif (yvarn=="tx") then iobs=3; imod=7 elseif (yvarn=="tm") then iobs=4; imod=8 endif it=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))==ireg) then Omax=Dat_mean(ilon,ilat,imon,iobs,ilead); Omin=Dat_std(ilon,ilat,imon,iobs,ilead) Fmax=Dat_mean(ilon,ilat,imon,imod,ilead); Fmin=Dat_std(ilon,ilat,imon,imod,ilead) Do iyr=1,Nyrs !-1 it=it+1 !X4d(ilon,ilat,1,iyr)=Y_inp(it) ! Obs !X4d(ilon,ilat,2,iyr)=Yest(it,1) ! Ann !X4d(ilon,ilat,3,iyr)=X_inp(it,1) ! Model if (Y_inp(it)/=Miss) X4d(ilon,ilat,1,iyr)= Y_inp(it)*(Omax-Omin) + Omin ! Obs if (Yest(it,1)/=Miss) X4d(ilon,ilat,2,iyr)= Yest(it,1)*(Fmax-Fmin) + Fmin ! Ann if (X4d(ilon,ilat,2,iyr)/=Miss.AND.X4d(ilon,ilat,2,iyr)<0) X4d(ilon,ilat,2,iyr)=0 if (X_inp(it,1)/=Miss) X4d(ilon,ilat,3,iyr)=X_inp(it,1)*(Fmax-Fmin) + Fmin ! Model Enddo endif Enddo Enddo !print*,"Size of Y_inp(:),X_inp(:,:): ",Size(Y_inp, Dim=1), Size(X_inp, Dim=1), Size(X_inp, Dim=2) !print*,"Convert: Nit=",it; pause End Subroutine Convert_to_lat_lon !-------------------------------------------------------------------------------- Subroutine Store_5th_Arr (ilead,imon) Integer :: imon,imem,ilon,ilat,it,iyr ! X6d(Nlon,Nlat,Nmon,Nyear,Nmem+2,4) - 1) OBS, 2) ANN est, 3) Model out, 4) Cli BC Do it=1,3 !iyr=1,Nyear Do iyr=1,Nyear ! imem=1,Nmem Do ilon=1,Nlon Do ilat=1,Nlat X5d(ilon,ilat,imon,iyr,it)=X4d(ilon,ilat,it,iyr) Enddo Enddo Enddo Enddo ! Climatological Bias Correction: Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))>0) then Do iyr=1,Nyear if (X4d(ilon,ilat,3,iyr)/=Miss.AND.Bias_Cli(ilon,ilat,ilead,imon)/=Miss) then if (yvarn=="pr") then X5d(ilon,ilat,imon,iyr,4)=X4d(ilon,ilat,3,iyr)*Bias_Cli(ilon,ilat,ilead,imon) ! For rainfall: BC = Model*BC_Factor else X5d(ilon,ilat,imon,iyr,4)=X4d(ilon,ilat,3,iyr)-Bias_Cli(ilon,ilat,ilead,imon) ! For Temp. X1d = Model-Bias = Corrected endif if (X5d(ilon,ilat,imon,iyr,4)<0) X5d(ilon,ilat,imon,iyr,4)=0 ! Monthly rainfall and Temp. must not be negative else X5d(ilon,ilat,imon,iyr,4)=Miss endif Enddo endif Enddo Enddo End Subroutine Store_5th_Arr !-------------------------------------------------------------------------------- Subroutine Make_Data (yvarn,ireg,imon,ilead) Character (Len=*) :: yvarn Integer :: ireg,imon,imon1,ilead,ipredictand,iprd Integer :: ilon,ilat,it,Npoint,Mx,i,j,Nspl,iyr Npoint = Get_Reg_Points (ireg) Ngrd = Npoint ! So diem luoi trong vung ireg Nsample1=Nyrs*Npoint ! Dung luong mau !print*,"Nsample1,Nyrs=",Nsample1,Nyrs Mx=18 ! So bien doc lap Mvars=Mx if (Allocated(X_inp)) Deallocate (X_inp) if (Allocated(Y_inp)) Deallocate (Y_inp) Allocate (X_inp(Nsample1,Mx),Y_inp(Nsample1)) if (yvarn=="pr") then ipredictand=1; iprd=5 elseif (yvarn=="t2") then ipredictand=2; iprd=6 elseif (yvarn=="tx") then ipredictand=3; iprd=7 elseif (yvarn=="tm") then ipredictand=4; iprd=8 endif i=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))==ireg) then Do it=1,Ntime imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon if (imon1==imon) then i=i+1 Y_inp(i) =Dat_All_0(ilon,ilat,it,ipredictand,ilead) ! Observation X_inp(i,1)=Dat_All_0(ilon,ilat,it,iprd,ilead) ! Model fcs X_inp(i,2)=Dat_All_0(ilon,ilat,it,11,ilead) ! olr X_inp(i,3)=Dat_All_0(ilon,ilat,it,12,ilead) ! q2 X_inp(i,4)=Dat_All_0(ilon,ilat,it,13,ilead) ! rh2 X_inp(i,5)=Dat_All_0(ilon,ilat,it,17,ilead) ! w700 !if (Dat_All(ilon,ilat,it,15,ilead)/=Miss.AND.Dat_All(ilon,ilat,it,16,ilead)/=Miss) then ! X_inp(i,6)=Dat_All(ilon,ilat,it,15,ilead)-Dat_All(ilon,ilat,it,16,ilead) ! u200-u850 !else ! X_inp(i,6)=Miss !endif X_inp(i,6)=Dat_All_0(ilon,ilat,it,16,ilead) ! u850 ! 4 neighbor gridpoints X_inp(i,7)=Dat_All_0(ilon-1,ilat,it,iprd,ilead) ! Model fcs X_inp(i,8)=Dat_All_0(ilon+1,ilat,it,iprd,ilead) ! Model fcs X_inp(i,9)=Dat_All_0(ilon,ilat-1,it,iprd,ilead) ! Model fcs X_inp(i,10)=Dat_All_0(ilon,ilat+1,it,iprd,ilead) ! Model fcs X_inp(i,11)=Dat_All_0(ilon-1,ilat,it,11,ilead) ! olr X_inp(i,12)=Dat_All_0(ilon+1,ilat,it,11,ilead) ! olr X_inp(i,13)=Dat_All_0(ilon,ilat-1,it,11,ilead) ! olr X_inp(i,14)=Dat_All_0(ilon,ilat+1,it,11,ilead) ! olr X_inp(i,15)=Dat_All_0(ilon-1,ilat,it,17,ilead) ! w700 X_inp(i,16)=Dat_All_0(ilon+1,ilat,it,17,ilead) ! w700 X_inp(i,17)=Dat_All_0(ilon,ilat-1,it,17,ilead) ! w700 X_inp(i,18)=Dat_All_0(ilon,ilat+1,it,17,ilead) ! w700 !if (Dat_All(ilon,ilat,it,9,ilead)/=Miss.AND.Dat_All(ilon,ilat,it,10,ilead)/=Miss) then ! X_inp(i,2)=Dat_All(ilon,ilat,it,9,ilead)-Dat_All(ilon,ilat,it,10,ilead) ! h500-h850 !else ! X_inp(i,2)=Miss !endif !X_inp(i,3)=Dat_All(ilon,ilat,it,11,ilead) ! olr !X_inp(i,4)=Dat_All(ilon,ilat,it,12,ilead) ! q2 !X_inp(i,5)=Dat_All(ilon,ilat,it,13,ilead) ! rh2 !X_inp(i,6)=Dat_All(ilon,ilat,it,14,ilead) ! slp !if (Dat_All(ilon,ilat,it,15,ilead)/=Miss.AND.Dat_All(ilon,ilat,it,16,ilead)/=Miss) then ! X_inp(i,7)=Dat_All(ilon,ilat,it,15,ilead)-Dat_All(ilon,ilat,it,16,ilead) ! u200-u850 !else ! X_inp(i,7)=Miss !endif !X_inp(i,8)=Dat_All(ilon,ilat,it,17,ilead) ! w700 endif Enddo endif Enddo Enddo Nspl=i 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 !print*,"Nsample1,i,Npoint,Nsample1/Npoint=",Nsample1,i,Npoint,Nsample1/Npoint; pause !Deallocate (Y_inp,X_inp) !Allocate (Y_inp(Nsample1),X_inp(Nsample1,Mx)) !Y_inp(1:Nsample1)=Outputs_Y(1,1:Nsample1); !Do j=1,Mx ! X_inp(1:Nsample1,j)=Inputs_X(j,1:Nsample1) !Enddo End Subroutine Make_Data !-------------------------------------------------------------------------------- Integer Function Get_Reg_Points (ireg) Integer :: ireg,ilon,ilat,Npoint Npoint=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)==ireg) Npoint=Npoint+1 Enddo Enddo Get_Reg_Points = Npoint End Function Get_Reg_Points !-------------------------------------------------------------------------------- Subroutine Calc_Bias_Cli (yvarn) Character (Len=*) :: yvarn !Integer, parameter :: Nt0=348, Nyr=Nt0/Nmon ! 1/1982-12/2010 Integer :: ilon,ilat,imon,imon1,it,ilead,iob,ifc,iyr !Real :: F(Nt0), O(Nt0), Xo(Nyr), Xf(Nyr), MeanO,MeanF Real :: F(Ntime),O(Ntime),Xo(Nyrs),Xf(Nyrs), MeanO,MeanF if (yvarn=="pr") then iob=1; ifc=5 elseif (yvarn=="t2") then iob=2; ifc=6 elseif (yvarn=="tx") then iob=3; ifc=7 elseif (yvarn=="tm") then iob=4; ifc=8 endif if (Allocated(Del_Cli)) Deallocate(Del_Cli) Allocate (Del_Cli(Nlon,Nlat,NLead+2,Nmon)); Del_Cli=Miss ! NLead+2: Mean Obs Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))>0) then Do ilead=1,NLead O(1:Ntime)=Dat_All(ilon,ilat,1:Ntime,iob,ilead) F(1:Ntime)=Dat_All(ilon,ilat,1:Ntime,ifc,ilead) Do imon=1,Nmon iyr=0; Do it=1,Ntime imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon if (imon1==imon) then iyr=iyr+1; Xo(iyr)=O(it); Xf(iyr)=F(it) endif Enddo MeanO=Ave(Xo,iyr,Miss); MeanF=Ave(Xf,iyr,Miss) if (MeanO/=Miss.AND.MeanF/=Miss.AND.MeanF/=0.) then if (yvarn=="pr") then ! For rainfall: Del_Cli=MeanO/MeanF Del_Cli(ilon,ilat,ilead,imon)=MeanO/MeanF ! Model_BC = Del_Cli*Model else Del_Cli(ilon,ilat,ilead,imon)=MeanF-MeanO ! Model_BC = MeanF - Del_Cli endif endif if (ilead==1) Del_Cli(ilon,ilat,NLead+2,imon)=MeanO ! NLead+2: Luu Mean Obs de su dung Enddo Enddo endif Enddo Enddo Del_Cli(:,:,NLead+1,:)=Del_Cli(:,:,NLead,:) ! Gan sai so cua LT=5 cho LT=6 End Subroutine Calc_Bias_Cli !-------------------------------------------------------------------------------- Subroutine Read_all_data (Dir,Model_nm) Character (Len=*) :: Dir,Model_nm Character (Len=255) :: Fname Real, Allocatable :: xlon(:),xlat(:),xtime(:) Integer :: ilead,ivar,Nlo,Nla,Nt Do ilead=1,NLead write(Fname,"(A,I1.1,A)") trim(Dir)//"lead",ilead,"_"//trim(Model_nm)//"_VN.nc" ! clWRF_VN.nc" Do ivar=1,Nvar Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,trim(vname(ivar))) if (ilead==1.AND.ivar==1) then Ntime=Nyrs*Nmon ! Ntime=Nt Allocate (Dat_All(Nlon,Nlat,Ntime,Nvar,NLead)) Dat_All=Miss endif Dat_All(:,:,1:Nt,ivar,ilead) = X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) Enddo Enddo End Subroutine Read_all_data !-------------------------------------------------------------------------------- Subroutine Data_standardization Integer :: ilon,ilat,imon,imon1,ilead,it,i,ivar, Nlead2 Real :: X1d(Nyrs),xmean,std,tmp, MaxX, MinX print*,"Data standardization..." Nlead2=Nlead+2 ! Lead=6,7 <-- Lead=5 Allocate (Dat_mean(Nlon,Nlat,Nmon,Nvar,NLead2),Dat_std(Nlon,Nlat,Nmon,Nvar,NLead2)) Dat_mean=Miss; Dat_std=Miss Allocate (Dat_All_0(Nlon,Nlat,Ntime,Nvar,NLead)) ! Stadardized data Dat_All_0=Miss Do ilon=1,Nlon Do ilat=1,Nlat Do ivar=1,Nvar Do ilead=1,Nlead Do imon=1,Nmon i=0 Do it=1,Ntime imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon if (imon1==imon) then i=i+1; X1d(i)=Dat_All(ilon,ilat,it,ivar,ilead) endif Enddo !X1d(1:Ntime)=Dat_All(ilon,ilat,1:Ntime,ivar,ilead) !xmean=Ave(X1d,i,Miss); std=Stdev(X1d,i,Miss) MaxX=Xmax(X1d,i,Miss); MinX=Xmin(X1d,i,Miss) Dat_mean(ilon,ilat,imon,ivar,ilead)=MaxX ! xmean Dat_std (ilon,ilat,imon,ivar,ilead)=MinX ! std Do it=1,Ntime imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon if (imon1==imon) then tmp = Dat_All(ilon,ilat,it,ivar,ilead) ! Unstandadized data !if (tmp/=Miss) then ! Dat_All_0(ilon,ilat,it,ivar,ilead)=(tmp-xmean)/std !endif if (tmp/=Miss.AND.MaxX/=Miss.AND.MinX/=Miss) then Dat_All_0(ilon,ilat,it,ivar,ilead)=(tmp-MinX)/(MaxX-MinX) endif endif Enddo Enddo Enddo Dat_mean(ilon,ilat,1:Nmon,ivar,Nlead+1)=Dat_mean(ilon,ilat,1:Nmon,ivar,Nlead) Dat_std (ilon,ilat,1:Nmon,ivar,Nlead+1)=Dat_std (ilon,ilat,1:Nmon,ivar,Nlead) Dat_mean(ilon,ilat,1:Nmon,ivar,Nlead+2)=Dat_mean(ilon,ilat,1:Nmon,ivar,Nlead) Dat_std (ilon,ilat,1:Nmon,ivar,Nlead+2)=Dat_std (ilon,ilat,1:Nmon,ivar,Nlead) Enddo Enddo Enddo End Subroutine Data_standardization !-------------------------------------------------------------------------------- !Include "/home/tanpv/Mylibs/Sub_HQTT.f90" !Include "/home/tanpv/Mylibs/Sub_LAD.f90" !Include "/home/tanpv/Mylibs/Sub_LMV.f90" Include "/home/tanpv/Mylibs/Calc_Correlation_Coefficient.f90" ! Correlation (X, Y, N, Corr, Prob, Miss) Include "/home/tanpv/Mylibs/Sub_Verification.f90" ! ME, MAE, RME, RMAE Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Sub_neural_2023.f90" Include "/home/tanpv/Mylibs/Stdev.f90" Include "/home/tanpv/Mylibs/Ave.f90" Include "/home/tanpv/Mylibs/Xmax.f90" Include "/home/tanpv/Mylibs/Xmin.f90" !-------------------------------------------------------------------------------- END Include "/home/tanpv/Mylibs/Student_Prob.f90" !-------------------------------------------------------------------------------------------