! ! 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=13, Nmem=6, Nfile=9 Integer, Parameter :: Nyrs=(Yr2-Yr1+1) Real, Parameter :: Miss=-99. Character (Len=4) :: vname(Nvar)=(/"pr","t2","tx","tm", & "h500","h850","olr","q2","rh2","slp","u200","u850","w700"/) ! 5 6 7 8 9 10 11 12 13 Integer :: Fcs_mon, Nyr1, NMax, Ncase, Ntc_var Character (Len=255) :: Fname Real, Allocatable :: lon(:),lat(:),Regs(:,:), X3d(:,:,:) Integer :: Nlon, Nlat, Ntime, Ngrd, ilon, ilat, it Integer :: Mvars, imem Real, Allocatable :: Dat_All(:,:,:,:,:) ! (Nlon,Nlat,NLead,Nvar,Nmem) - Anomaly data = Original data - Original Clim Real, Allocatable :: Dat_Org(:,:,:,:,:) ! (Nlon,Nlat,NLead,Nvar,Nmem) - Original data Real, Allocatable :: Clim (:,:,:,:,:) ! (Nlon,Nlat,Nmon,Nvar,NLead) - Original Clim 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 !-------------------------------------------------------------------------------- ! Regression: Real, Allocatable :: Xorg(:), Xhq(:) ! Real, Allocatable :: Beta(:,:,:,:) ! NLead,Nmon,Nreg,Mhq Integer :: Mhq, Nhq ! Mhq = Mxvar+1, Nhq=Nsample*Mhq 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, Mask_file, Dir_ANN_Out, Dir_Cli, Dir_ANN_Save, Dir_BC_Ope, Dir_Corrected 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 !-------------------------------------------------------------------------------- Integer :: Op_Nhidden, Op_Nepoc ! Optimal values Real :: Op_learning_rate Real, Allocatable :: X4d(:,:,:,:), F1d(:),O1d(:), X5d(:,:,:,:),Score_Cor(:,:,:), Score_AV(:,:,:) Real, Allocatable :: xlon(:), xlat(:), xtime(:), xlev(:) Real, Allocatable :: X6d(:,:,:,:,:), BC_Cli(:,:,:,:) Integer :: Nlo, Nla, Nt, Nz, yyyy, mm !------------------- For ANN ------------------------------------------------------------- Namelist /ann_info/ Nhidden,learning_rate,momentum,training_mode,Nepoc,debug Namelist /Dir_info/ Mask_file,Dir_ANN_Save,Dir_Cli,Dir_Corrected,Dir_BC_Ope Namelist /Model_info/ Model_nm, yvarn Namelist /Time_info/ yyyy, mm ! Nhidden = 5 !4 ! Number of hidden layers ! learning_rate = 0.01 ! momentum = 0.9 ! 10.0 ! training_mode = .false. ! For training ! Nepoc = 1000 ! Iterative training times ! debug = .false. !-------------------------------------------------------------------------------- Open (1, file="Namelist_BC_Ope.txt",Status="old") Read (1, ann_info) Read (1, Dir_info) Read (1, Model_info) Read (1, Time_info) learning_rate0 = learning_rate training_mode = .false. !Fname="/home/tanpv/Data/AWC/Mask_VN_0.25_Corrected.nc" !Fname="/home/tanpv/Progs/Dr_VN_Ope/DATA/Mask_VN_0.25_Corrected.nc" Call rd_ncdf_2d (trim(Mask_file),lon,lat,Nlon,Nlat,Regs,"reg") ! 7 sub-regs mask !Dir_ANN_Save ="/home/tanpv/Data/562_VNU/0.25_VN/2024/ANN_Save/" !Dir_Cli ="/work/projects/G09/562_VNU/0.25_VN/2024/Cli_Out/" !"/home/tanpv/Data/562_VNU/0.25_VN/Cli_Out/" !Dir_Corrected ="/home/tanpv/Data/562_VNU/0.25_VN/2024/Corrected_Ope/" !Dir_ANN_Save ="/home/tanpv/Progs/Dr_VN_Ope/DATA/ANN_Save/" !Dir_Cli ="/home/tanpv/Progs/Dr_VN_Ope/DATA/Cli_Out/" !"/home/tanpv/Data/562_VNU/0.25_VN/Cli_Out/" !Dir_Corrected ="/home/tanpv/Progs/Dr_VN_Ope/DATA/Corrected_Ope/" !-------------------------------------------------------------------------------- write (Fname,"(A)") trim(Dir_Cli)//"Clim_Mean_Bias_"//trim(Model_nm)//"_"//trim(yvarn)//".nc" Call Read_Clim_Bias (Fname) ! Climatological Mean and Biases Fname=trim(Dir_Cli)//"Beta_"//trim(yvarn)//"_"//trim(Model_nm)//".nc" Call Read_Beta_HQ (Fname) ! Regression Coefficiences !Write (Dir,"(A,I4.4,A,I2.2,A)") "/home/tanpv/Data/562_VNU/0.25_VN/For_BC_Ope/",yyyy,"/",mm,"/" !Write (Dir,"(A,I4.4,A,I2.2,A)") "/home/tanpv/Progs/Dr_VN_Ope/DATA/For_BC_Ope/",yyyy,"/",mm,"/" Write (Dir,"(A,I4.4,A,I2.2,A)") trim(Dir_BC_Ope),yyyy,"/",mm,"/" Call Read_all_data (Dir,Model_nm) !,imem) Call Calc_Anomaly Nyear=Nyrs if (Allocated(X6d)) Deallocate (X6d) Allocate (X6d(Nlon,Nlat,NLead,Nmem+1,4)) ! 4: 1) ANN, 2) Mod, 3) Cli, 4) MLR X6d=Miss Do ilead=1,NLead imon=mm+ilead ! To get data for ANN and MLR if (imon>Nmon) imon=imon-Nmon ! imon of the year if (Allocated(X5d)) Deallocate (X5d); Allocate (X5d(Nlon,Nlat,3,Nmem)) ! In X5d z=3: 1) ann, 2) model, 3) mlr X5d=Miss Do ireg=1,Nreg print*,"ireg,imon,ilead=",ireg,imon,ilead Do imem=1,Nmem Call Make_Data (trim(yvarn),ireg,imon,ilead,imem) 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 print*,"Calling ANN..." !Call ANN_2023(Inputs_X,Outputs_Y,Nsample1,Mvars,1,Miss) Call ANN_2024(Inputs_X,Outputs_Y,Nsample1,Mvars,1,Miss) Call Calc_HQTT (ilead,imon,ireg) ! Return to Xhq(Nsample1) Call Convert_to_lat_lon (ireg,imem) ! Store in X5d Enddo Enddo Call Store_6th_Arr (ilead,imon) Enddo write (Fname,"(A,I4.4,I2.2,A)") trim(Dir_Corrected)//"BC_Ope_",yyyy,mm,"_"//trim(yvarn)//"_"//trim(Model_nm)//".nc" !"_lead_",ilead,".nc" Call Write_Out_4_Drought_Ope (Fname,yvarn) !-------------------------------------------------------------------------------- CONTAINS Subroutine Calc_HQTT (ilead,imon,ireg) Integer :: i,j,k, ilead,imon,ireg Real :: S,Q,U,RR,F, tmp Real :: xtmp(Mvars) Logical :: OK ! Tinh gia tri HQ: Beta(ilead,imon,ireg,1:Mhq) if (Allocated (Xhq)) Deallocate (Xhq); Allocate (Xhq(Nsample1)); Xhq=Miss Do i=1,Nsample1 xtmp(1:Mvars)=X_inp(i,1:Mvars); OK=Valid_Ar1(xtmp,Mvars,Miss) !if (Y_inp(i)/=Miss.AND.OK) then if (OK) then tmp=0. Do j=2,Mhq tmp=tmp+Beta(ilead,imon,ireg,j) * X_inp(i,j-1) Enddo Xhq(i)=Beta(ilead,imon,ireg,1) + tmp else Xhq(i)=Miss endif Enddo End Subroutine Calc_HQTT !-------------------------------------------------------------------------------- Subroutine Write_Out_4_Drought_Ope (Fname,yvarn) Character (Len=*) :: Fname,yvarn Integer, parameter :: Nvar=4 ! for wrting out 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)=(/"ann","mod","cli","mlr"/) Character (Len=20), parameter :: v_ln(Nvar)=(/"ANN estimation","Model Forecast","Climatology bias correction","MLR estimation"/) ! X6d(Nlon,Nlat,NLead,Nmem+1,4) - 1) ANN, 2) Mod, 3) Cli, 4) MLR Nz=Nmem+1; Allocate (xtime(NLead),zlev(Nz),Y4d(Nlon,Nlat,Nz,NLead)) Do i=1,Nz; zlev(i)=i; Enddo Do i=1,NLead; xtime(i)=i; Enddo write (t_unit,"(A,I4.4,A,I2.2,A)") "months since ",yyyy,"-",mm+1,"-01" lon_nm="lon"; lat_nm="lat"; z_nm="member"; t_nm="lead" var_ln="Outputs of "//trim(yvarn)//" for drought calculation. t=1..6: Member 1..6; t=7: Ens. " Do ivar=1,Nvar Do imem=1,Nz Do ilead=1,NLead Y4d(:,:,imem,ilead)=X6d(:,:,ilead,imem,ivar) Enddo Enddo var_longname=trim(var_ln)//v_ln(ivar) if (ivar==1) then Call wr_ncdf_4d (Fname,Lon,Lat,zlev,xtime,Nlon,Nlat,Nz,NLead,Y4d,trim(vname(ivar)), & t_unit, var_longname, "-", lon_nm, lat_nm, z_nm, t_nm, Miss) else Call wr_ncdf_4d_apd (Fname,Y4d,trim(vname(ivar)),Nlon,Nlat,Nz,NLead, & var_longname,"-", lon_nm, lat_nm, z_nm, t_nm, Miss) endif Enddo End Subroutine Write_Out_4_Drought_Ope !-------------------------------------------------------------------------------- Subroutine Store_6th_Arr (ilead,imon) Integer :: ilead,imon,imem,ilon,ilat,it,iyr,iob,ifc Real :: o1,o2,o3,o4, c1,c2,c3,c4, bias Real :: X1d(Nmem) ! X6d(Nlon,Nlat,NLead,Nmem+1,4) - 1) ANN, 2) Mod, 3) Cli, 4) MLR ! X5d(Nlon,Nlat,3,Nmem) 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 Do imem=1,Nmem Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))>0) then o1=X5d(ilon,ilat,1,imem); ! Ann Ano o2=X5d(ilon,ilat,2,imem); ! Mod Ano o3=X5d(ilon,ilat,3,imem) ! Mlr Ano c1=Clim(ilon,ilat,imon, iob,ilead); ! Tham so imon chi de xac dinh trung binh khi hau c2=Clim(ilon,ilat,imon, ifc,ilead); c3=c1 if (o1/=Miss.AND.c1/=Miss) X6d(ilon,ilat,ilead,imem,1) = o1+c1 ! Ann Ano + Obs Clim = Ann if (o2/=Miss.AND.c2/=Miss) X6d(ilon,ilat,ilead,imem,2) = o2+c2 ! Mod Ano + Mod Clim = Org Mod if (o3/=Miss.AND.c3/=Miss) X6d(ilon,ilat,ilead,imem,4) = o3+c3 ! Mlr Ano + Obs Clim = MLR <<<< Chu y chi so 4 ! Bias Cli correction: Bias_Cli(ilon,ilat,ilead,imon) bias=Bias_Cli(ilon,ilat,ilead,imon) if (o2/=Miss.AND.bias/=Miss) then if (yvarn=="pr") then X6d(ilon,ilat,ilead,imem,3) = X6d(ilon,ilat,ilead,imem,2) * bias ! Model correction else X6d(ilon,ilat,ilead,imem,3) = X6d(ilon,ilat,ilead,imem,2) + bias endif if (X6d(ilon,ilat,ilead,imem,3)<0) & X6d(ilon,ilat,ilead,imem,3) = X6d(ilon,ilat,ilead,imem,2) ! Must not be negative, assign to model value else X6d(ilon,ilat,ilead,imem,3) = Miss endif if (X6d(ilon,ilat,ilead,imem,1)<0) & X6d(ilon,ilat,ilead,imem,1) = X6d(ilon,ilat,ilead,imem,3) ! Must not be negative, assign to Cli value if (X6d(ilon,ilat,ilead,imem,4)<0) & X6d(ilon,ilat,ilead,imem,4) = X6d(ilon,ilat,ilead,imem,3) ! Must not be negative, assign to Cli value endif Enddo Enddo Enddo !print*,X6d(:,:,imon,1,1); pause !Enddo ! Do ensemble and assign to the position of Nmem+1: !Do iyr=1,Nyear Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))>0) then X1d(1:Nmem)=X6d(ilon,ilat,ilead,1:Nmem,1); X6d(ilon,ilat,ilead,Nmem+1,1)=Ave(X1d,Nmem,Miss) X1d(1:Nmem)=X6d(ilon,ilat,ilead,1:Nmem,2); X6d(ilon,ilat,ilead,Nmem+1,2)=Ave(X1d,Nmem,Miss) X1d(1:Nmem)=X6d(ilon,ilat,ilead,1:Nmem,3); X6d(ilon,ilat,ilead,Nmem+1,3)=Ave(X1d,Nmem,Miss) X1d(1:Nmem)=X6d(ilon,ilat,ilead,1:Nmem,4); X6d(ilon,ilat,ilead,Nmem+1,4)=Ave(X1d,Nmem,Miss) endif Enddo Enddo !Enddo End Subroutine Store_6th_Arr !-------------------------------------------------------------------------------- Subroutine Convert_to_lat_lon (ireg,imem) Integer :: ireg,imem,ilon,ilat,it,iyr it=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))==ireg) then !Do iyr=1,Nyear it=it+1 !X5d(ilon,ilat,1,iyr,imem)=Y_inp(it) ! OBS Anomaly X5d(ilon,ilat,1,imem)=Yest(it,1) ! ANN Anomaly estimation X5d(ilon,ilat,2,imem)=X_inp(it,1) ! Origin model Anomaly X5d(ilon,ilat,3,imem)=Xhq(it) ! MLR Anomaly estimation !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*,"Nsample1,Convert: Nit=",Nsample1,it; pause End Subroutine Convert_to_lat_lon !-------------------------------------------------------------------------------- Subroutine Make_Data (yvarn,ireg,imon,ilead,imem) Character (Len=*) :: yvarn Integer :: ireg,imon,imon1,ilead,ipredictand,iprd Integer :: ilon,ilat,it,Npoint,Mx,i,j,Nspl,imem Npoint = Get_Reg_Points (ireg) Ngrd = Npoint ! So diem luoi trong vung ireg Nsample1=Nyrs*Npoint ! Dung luong mau !print*,"Npoint,Nsample1,Nyrs,Ntime=",Npoint,Nsample1,Nyrs,Ntime print*,"Beginning Make_Data" Mx=22 ! 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 ! ilead tuong duong voi Initial montn (mm) + 1 i=i+1 Y_inp(i) =Miss !Dat_All(ilon,ilat,it,ipredictand,ilead,imem) ! Observation X_inp(i,1)=Dat_All(ilon,ilat,ilead,iprd,imem) ! Model fcs X_inp(i,2)=Dat_All(ilon,ilat,ilead,7,imem) ! olr X_inp(i,3)=Dat_All(ilon,ilat,ilead,8,imem) ! q2 X_inp(i,4)=Dat_All(ilon,ilat,ilead,9,imem) ! rh2 X_inp(i,5)=Dat_All(ilon,ilat,ilead,13,imem) ! w700 if (Dat_All(ilon,ilat,ilead,11,imem)/=Miss.AND.Dat_All(ilon,ilat,ilead,12,imem)/=Miss) then X_inp(i,6)=Dat_All(ilon,ilat,ilead,11,imem)-Dat_All(ilon,ilat,ilead,12,imem) ! u200-u850 else X_inp(i,6)=Miss endif ! 4 neighbor gridpoints X_inp(i,7)=Dat_All(ilon-1,ilat,ilead,iprd,imem) ! Model fcs X_inp(i,8)=Dat_All(ilon+1,ilat,ilead,iprd,imem) ! Model fcs X_inp(i,9)=Dat_All(ilon,ilat-1,ilead,iprd,imem) ! Model fcs X_inp(i,10)=Dat_All(ilon,ilat+1,ilead,iprd,imem) ! Model fcs X_inp(i,11)=Dat_All(ilon-1,ilat,ilead,7,imem) ! olr X_inp(i,12)=Dat_All(ilon+1,ilat,ilead,7,imem) ! olr X_inp(i,13)=Dat_All(ilon,ilat-1,ilead,7,imem) ! olr X_inp(i,14)=Dat_All(ilon,ilat+1,ilead,7,imem) ! olr X_inp(i,15)=Dat_All(ilon-1,ilat,ilead,13,imem) ! w700 X_inp(i,16)=Dat_All(ilon+1,ilat,ilead,13,imem) ! w700 X_inp(i,17)=Dat_All(ilon,ilat-1,ilead,13,imem) ! w700 X_inp(i,18)=Dat_All(ilon,ilat+1,ilead,13,imem) ! w700 X_inp(i,19)=Dat_All(ilon-1,ilat,ilead,8,imem) ! q2 X_inp(i,20)=Dat_All(ilon+1,ilat,ilead,8,imem) ! q2 X_inp(i,21)=Dat_All(ilon,ilat-1,ilead,8,imem) ! q2 X_inp(i,22)=Dat_All(ilon,ilat+1,ilead,8,imem) ! q2 !if (Dat_All(ilon,ilat,it,9,ilead,imem)/=Miss.AND.Dat_All(ilon,ilat,it,10,ilead,imem)/=Miss) then ! X_inp(i,2)=Dat_All(ilon,ilat,it,9,ilead,imem)-Dat_All(ilon,ilat,it,10,ilead,imem) ! h500-h850 !else ! X_inp(i,2)=Miss !endif !X_inp(i,3)=Dat_All(ilon,ilat,it,11,ilead,imem) ! olr !X_inp(i,4)=Dat_All(ilon,ilat,it,12,ilead,imem) ! q2 !X_inp(i,5)=Dat_All(ilon,ilat,it,13,ilead,imem) ! rh2 !X_inp(i,6)=Dat_All(ilon,ilat,it,14,ilead,imem) ! slp !if (Dat_All(ilon,ilat,it,15,ilead,imem)/=Miss.AND.Dat_All(ilon,ilat,it,16,ilead,imem)/=Miss) then ! X_inp(i,7)=Dat_All(ilon,ilat,it,15,ilead,imem)-Dat_All(ilon,ilat,it,16,ilead,imem) ! u200-u850 !else ! X_inp(i,7)=Miss !endif !X_inp(i,8)=Dat_All(ilon,ilat,it,17,ilead,imem) ! w700 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"(A,8I5)","Nspl,Nsample1,i,Npoint,Nsample1/Npoint=",Nspl,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 Do i=1,Nsample1 X_inp(i,j)=Inputs_X(j,i) Enddo Enddo print*,"End of Make_Data" 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_Anomaly Integer :: ilon,ilat,imon,imon1,it,ilead,iob,ifc,iyr,ireg,imem Real :: Mean, xtmp ! Dat_Org(Nlon,Nlat,NLead,Nvar,Nmem) Allocate (Dat_All(Nlon,Nlat,NLead,Nvar,Nmem)); Dat_All=Miss ! Anomaly Do ilead=1,NLead imon=mm+ilead ! imon = Initial month (mm) + ilead if (imon>Nmon) imon=imon-Nmon !print*,imon,ilead Do imem=1,Nmem Do ivar=1,Nvar Do ilon=1,Nlon ! All domain Do ilat=1,Nlat Mean = Clim(ilon,ilat,imon,ivar+4,ilead) ! ivar=1-4: Obs (pr_o, t2_o, tx_o, tm_o) xtmp = Dat_Org(ilon,ilat,ilead,ivar,imem) if (Mean/=Miss.AND.xtmp/=Miss) Dat_All(ilon,ilat,ilead,ivar,imem)=xtmp-Mean Enddo Enddo Enddo Enddo Enddo End Subroutine Calc_Anomaly !-------------------------------------------------------------------------------- Subroutine Read_all_data (Dir,Model_nm) Character (Len=*) :: Dir,Model_nm Character (Len=4), parameter :: varn(Nvar)=(/"pr","tas","tasmax","tasmin", & "h500","h850","olr","q2","rh2","slp","u200","u850","w700"/) Character (Len=255) :: Fname Real, Allocatable :: xlon(:),xlat(:),xtime(:) Integer :: imem,ilead,ivar,Nlo,Nla,Nt, k0, it1, it2 Do imem=1,Nmem write(Fname,"(A,I2.2,A)") trim(Dir)//trim(Model_nm)//"_",imem,".nc" Do ivar=1,Nvar Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,trim(varn(ivar))) if (imem==1.AND.ivar==1) then Allocate (Dat_Org(Nlon,Nlat,NLead,Nvar,Nmem)); Dat_Org=Miss endif Dat_Org(:,:,:,ivar,imem) = X3d(:,:,:) Deallocate (xlon,xlat,xtime,X3d) Enddo Enddo End Subroutine Read_all_data !-------------------------------------------------------------------------------- Subroutine Read_Beta_HQ (Fname) Character (Len=*) :: Fname Integer :: imem,ilead,ivar,Nlo,Nla,Nz,Nt,imon, k0, it1, it2 Real, Allocatable :: xlon(:),xlat(:),xlev(:),xtime(:) Call rd_ncdf_4d (Fname,xlon,xlat,xlev,xtime,Nlo,Nla,Nz,Nt,Beta,"beta") Mhq=Nt Deallocate (xlon,xlat,xlev,xtime) End Subroutine Read_Beta_HQ !-------------------------------------------------------------------------------- Subroutine Read_Clim_Bias (Fname) Character (Len=*) :: Fname Integer, Parameter :: Nvar4=Nvar+4 Real, Allocatable :: xlon(:),xlat(:),xlev(:),xtime(:),X4d(:,:,:,:) Integer :: imem,ilead,ivar,Nlo,Nla,Nz,Nt,imon, k0, it1, it2 Character (Len=4) :: vname(Nvar4)=(/"pr_o","t2_o","tx_o","tm_o","pr","t2","tx","tm", & "h500","h850","olr","q2","rh2","slp","u200","u850","w700"/) ! Clim (:,:,:,:,:) ! (Nlon,Nlat,Nmon,Nvar,NLead) !write (Fname,"(A)") trim(Dir_Cli)//"Clim_Mean_Bias_"//trim(Model_nm)//"_"//trim(yvarn)//".nc" Allocate (Clim(Nlon,Nlat,Nmon,Nvar4,NLead)) Call rd_ncdf_4d (Fname,xlon,xlat,xlev,xtime,Nlo,Nla,Nz,Nt,Bias_Cli,"bias") ! Reading Bias_Cli = Model_Fcs - Obs (For rainfall: Bias_Cli=Obs/Model_Fcs Deallocate (xlon,xlat,xlev,xtime) Do ivar=1,Nvar4 Call rd_ncdf_4d (Fname,xlon,xlat,xlev,xtime,Nlo,Nla,Nz,Nt,X4d,trim(vname(ivar))) !print*,trim(vname(ivar)) Do imon=1,Nmon Do ilead=1,NLead Clim(:,:,imon,ivar,ilead)=X4d(:,:,ilead,imon) Enddo Enddo Deallocate (xlon,xlat,xlev,xtime,X4d) Enddo End Subroutine Read_Clim_Bias !-------------------------------------------------------------------------------- Logical Function Miss_Ar1(X,N,Miss) Integer :: N,i Real :: X(N), Miss Logical :: Empty Empty=.False. Do i=1,N if (X(i)==Miss) then Empty=.True.; Exit endif Enddo Miss_Ar1=Empty End Function Miss_Ar1 !-------------------------------------------------------------------------------- !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/Sub_neural_2024.f90" Include "/home/tanpv/Mylibs/Stdev.f90" Include "/home/tanpv/Mylibs/Ave.f90" Include "/home/tanpv/Mylibs/Function_Valid_Ar1.f90" !-------------------------------------------------------------------------------- END Include "/home/tanpv/Mylibs/Student_Prob.f90" !-------------------------------------------------------------------------------------------