!
!  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=17, Nmem=6, 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(:,:,:)
   Integer             :: Nlon, Nlat, Ntime, Ngrd, ilon, ilat, it
   Integer             :: Mvars, imem
   Real, Allocatable   :: Dat_All  (:,:,:,:,:,:) ! (Nlon,Nlat,Ntime,Nvar,NLead,Nmem)
   Real, Allocatable   :: Dat_All_0(:,:,:,:,:,:) ! (Nlon,Nlat,Ntime,Nvar,NLead,Nmem)
   Real, Allocatable   :: Dat_mean(:,:,:,:,:)  ! Dat_mean(Nlon,Nlat,Nmon,Nvar,NLead),Dat_std(Nlon,Nlat,Nmon,Nvar,NLead)
   Real, Allocatable   :: Dat_std (:,:,:,:,:)
   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)
   !Real, Allocatable, Dimension(:,:,:,:) :: Pr_CORbc, Pr_Pvalbc, Pr_AVbc   ! between ANN and BC

   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
   !--------------------------------------------------------------------------------
   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
   !------------------- 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 = .false.       ! 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
   training_mode  = .false.

   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"
!   Do imem=1,1 !Nmem

   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/"
   !--------------------------------------------------------------------------------

      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")  ! Reading Bias_Cli = Model_Fcs - Obs (For rainfall: Bias_Cli=Obs/Model_Fcs

      Dir="/home/tanpv/Progs/562_VNU/Data/For_BC/"
      Call Read_all_data (Dir,Model_nm) !,imem)
      Call Data_standardization
 
      Nyear=Nyrs
      Do ilead=1,NLead
         if (Allocated(X6d)) Deallocate (X6d)
         Allocate (X6d(Nlon,Nlat,Ntime,Nmem+1,4))   ! Nmem+1: Ens; 6th-dimension = 4: 1) OBS, 2) ANN est, 3) Model out 4) Cli BC 
         !Do imon=1,Nmon
            if (Allocated(X4d)) Deallocate (X4d)
            if (Allocated(X5d)) Deallocate (X5d); 
            if (Allocated(BC_Cli)) Deallocate(BC_Cli)
            Allocate (BC_Cli(Nlon,Nlat,Ntime,Nmem)); BC_Cli=Miss
            Allocate (X4d(Nlon,Nlat,4,Ntime),X5d(Nlon,Nlat,3,Ntime,Nmem)) ! Nyrs-1))  ! z=4: 1) Yobs, 2) Yest, 3) Ymodel, 4) Ymodel-Bias
            X4d=Miss; 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)
                  if (ilead<NLead) then
                     !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"
                     write (save_file,"(2(A,I2.2),A)") trim(Dir_ANN_Save)//"Sav_"//trim(yvarn)//"_"//trim(Model_nm)//"_Reg_",ireg,"_lead_",ilead,".txt"
                  elseif (ilead==NLead) then
                     !write (save_file,"(3(A,I2.2),A)") trim(Dir_ANN_Save)//"Sav_"//trim(yvarn)//"_"//trim(Model_nm)//"_Reg_",ireg,"_mon_",imon,"_lead_",ilead-1,".txt"
                     write (save_file,"(2(A,I2.2),A)") trim(Dir_ANN_Save)//"Sav_"//trim(yvarn)//"_"//trim(Model_nm)//"_Reg_",ireg,"_lead_",ilead-1,".txt"
                  endif
                  learning_rate = learning_rate0  !  0.05
                  print*,"Calling ANN..."
                  Call ANN_2023(Inputs_X,Outputs_Y,Nsample1,Mvars,1,Miss)
                  Call Convert_to_lat_lon (ireg,ilead,imem)
               Enddo
            Enddo
            Call Ens_all_members (ilead) !,imon)
            !Call Verification (ilead,imon)
            !Call Write_Out (output_file,yvarn)
            Call Store_6th_Arr !(imon)
         !Enddo
         write (Fname,"(A,I2.2,A)") trim(Dir_4Dr)//"Dr_testing_"//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=Nmem+1; 
      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 2011-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..6: Member 1..6; Z=7: Ens. "
      Do ivar=1,Nvar
         Do imem=1,Nz
            !it=0
            !Do iyr=1,Nyear
            !   Do imon=1,Nmon
            !      it=it+1; Y4d(:,:,imem,it)=X6d(:,:,imon,iyr,imem,ivar)
            !   Enddo
            !Enddo
            Do it=1,Ntime
               Y4d(:,:,imem,it)=X6d(:,:,it,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,Ntime,Y4d,vname(ivar), &
                       t_unit, var_longname, "-", lon_nm, lat_nm, z_nm, t_nm, Miss)
         else
            Call wr_ncdf_4d_apd (Fname,Y4d,vname(ivar),Nlon,Nlat,Nz,Ntime, &
                           var_longname,"-", lon_nm, lat_nm, z_nm, t_nm, Miss)
         endif
      Enddo
   End Subroutine Write_Out_4_Drought
   !--------------------------------------------------------------------------------
   Subroutine Store_6th_Arr ! (imon)
      Integer :: imon,imem,ilon,ilat,it,iyr,ivar
      ! X6d(Nlon,Nlat,Nmon,Nyear,Nmem+2,4) - 1) OBS, 2) ANN est, 3) Model out, 4) Cli BC
      Do ivar=1,3             ! 1) OBS, 2) ANN est, 3) Model out
         Do iyr=1,Ntime !Nyear    
            Do imem=1,Nmem
               X6d(:,:,iyr,imem,ivar)=X5d(:,:,ivar,iyr,imem)
            Enddo
            X6d(:,:,iyr,Nmem+1,ivar)=X4d(:,:,ivar,iyr)
         Enddo
      Enddo
      ! it=4: BC Cli
      Do iyr=1,Ntime ! Nyear
         Do imem=1,Nmem
            X6d(:,:,iyr,imem,4)=BC_Cli(:,:,iyr,imem)  ! BC Cli members
         Enddo
         X6d(:,:,iyr,Nmem+1,4)=X4d(:,:,4,iyr)         ! BC Cli ensemble mean
      Enddo
   End Subroutine Store_6th_Arr
   !--------------------------------------------------------------------------------
   Subroutine Ens_all_members (ilead) !,imon)
      Integer :: imon,imon1,ireg,imem,ilon,ilat,it,iyr
      Real :: X1d(Nmem)

      Do ilon=1,Nlon
         Do ilat=1,Nlat
            if (Int(Regs(ilon,ilat))>0) then
               Do iyr=1,Ntime !Nyear ! 7 ! Nyrs-1 ==============> SL moi chi den 2018
                  imon=mod(iyr,Nmon); if (imon==0) imon=Nmon
                  X1d(1:Nmem)=X5d(ilon,ilat,1,iyr,1:Nmem)    ! OBS 
                  X4d(ilon,ilat,1,iyr)=Ave(X1d, Nmem, Miss)
                  X1d(1:Nmem)=X5d(ilon,ilat,2,iyr,1:Nmem)    ! ANN out
                  X4d(ilon,ilat,2,iyr)=Ave(X1d, Nmem, Miss)
                  X1d(1:Nmem)=X5d(ilon,ilat,3,iyr,1:Nmem)    ! Model out
                  X4d(ilon,ilat,3,iyr)=Ave(X1d, Nmem, Miss)
                  Do imem=1,Nmem
                     !Do imon=1,Nmon
                        !imon=mod(iyr,Nmon); if (imon==0) imon=Nmon
                        !if (imon1==imon) then
                           if (X5d(ilon,ilat,3,iyr,imem)/=Miss.AND.Bias_Cli(ilon,ilat,ilead,imon)/=Miss) then
                              if (yvarn=="pr") then
                                 X1d(imem)=X5d(ilon,ilat,3,iyr,imem)*Bias_Cli(ilon,ilat,ilead,imon)      ! For rainfall: BC_Cli=Model*Bias
                              else
                                 X1d(imem)=X5d(ilon,ilat,3,iyr,imem)-Bias_Cli(ilon,ilat,ilead,imon)      ! for temp. BC_Cli = Model-Bias = Corrected
                              endif
                              if (X1d(imem)<0) X1d(imem)=0          ! Monthly values must not be negative
                              BC_Cli(ilon,ilat,iyr,imem)=X1d(imem)  ! BC Cli for each member
                           else
                              X1d(imem)=Miss; BC_Cli(ilon,ilat,iyr,imem)=Miss
                           endif
                        !endif
                     !Enddo
                  Enddo
                  X4d(ilon,ilat,4,iyr)=Ave(X1d, Nmem, Miss)  ! Ensemble mean of BC Cli
               Enddo
            endif
         Enddo
      Enddo
      ! Bias_Cli(:,:,:,:,:)  ! (Nlon,Nlat,NLead,Nmon)
   End Subroutine Ens_all_members
   !--------------------------------------------------------------------------------
   Subroutine Convert_to_lat_lon (ireg,ilead,imem)
      Integer :: imon,ireg,imem,ilon,ilat,it,iyr
      Real    :: MeanO, MeanF, stdO, stdF

      ! Dat_mean(Nlon,Nlat,Nvar,NLead),Dat_std(Nlon,Nlat,Nvar,NLead) xa=(x-mean)/std
      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
               !MeanO=Dat_mean(ilon,ilat,iobs,ilead); stdO=Dat_std(ilon,ilat,iobs,ilead)
               !MeanF=Dat_mean(ilon,ilat,imod,ilead); stdF=Dat_std(ilon,ilat,imod,ilead)
               Do iyr=1,Ntime ! Nyear
                  imon=mod(iyr,Nmon); if (imon==0) imon=Nmon
                  MeanO=Dat_mean(ilon,ilat,imon,iobs,ilead); stdO=Dat_std(ilon,ilat,imon,iobs,ilead)
                  MeanF=Dat_mean(ilon,ilat,imon,imod,ilead); stdF=Dat_std(ilon,ilat,imon,imod,ilead)
                  it=it+1
                  !if (Y_inp(it)/=Miss)   X5d(ilon,ilat,1,iyr,imem)=  Y_inp(it)*stdO + MeanO     ! OBS
                  !if (Yest(it,1)/=Miss)  X5d(ilon,ilat,2,iyr,imem)= Yest(it,1)*stdF + MeanF     ! ANN estimation
                  !if (X5d(ilon,ilat,2,iyr,imem)/=Miss.AND.X5d(ilon,ilat,2,iyr,imem)<0) X5d(ilon,ilat,2,iyr,imem)=0
                  !if (X_inp(it,1)/=Miss) X5d(ilon,ilat,3,iyr,imem)=X_inp(it,1)*stdF + MeanF     ! Origin model fcs

                  if (Y_inp(it)/=Miss)   X5d(ilon,ilat,1,iyr,imem)=  Y_inp(it)*(MeanO-stdO) + StdO     ! OBS
                  if (Yest(it,1)/=Miss)  X5d(ilon,ilat,2,iyr,imem)= Yest(it,1)*(MeanF-stdF) + StdF     ! ANN estimation
                  if (X5d(ilon,ilat,2,iyr,imem)/=Miss.AND.X5d(ilon,ilat,2,iyr,imem)<0) X5d(ilon,ilat,2,iyr,imem)=0
                  if (X_inp(it,1)/=Miss) X5d(ilon,ilat,3,iyr,imem)=X_inp(it,1)*(MeanF-stdF) + StdF     ! Origin model fcs
               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 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
      Nsample1=Ntime*Npoint

      !print*,"Npoint,Nsample1,Nyrs,Ntime=",Npoint,Nsample1,Nyrs,Ntime
      print*,"Beginning Make_Data"
      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,imem)   ! Observation
                     X_inp(i,1)=Dat_All_0(ilon,ilat,it,iprd,ilead,imem)          ! Model fcs
                     X_inp(i,2)=Dat_All_0(ilon,ilat,it,11,ilead,imem)   ! olr
                     X_inp(i,3)=Dat_All_0(ilon,ilat,it,12,ilead,imem)   ! q2
                     X_inp(i,4)=Dat_All_0(ilon,ilat,it,13,ilead,imem)   ! rh2
                     X_inp(i,5)=Dat_All_0(ilon,ilat,it,17,ilead,imem)   ! w700
                     !if (Dat_All_0(ilon,ilat,it,15,ilead,imem)/=Miss.AND.Dat_All_0(ilon,ilat,it,16,ilead,imem)/=Miss) then
                     !   X_inp(i,6)=Dat_All_0(ilon,ilat,it,15,ilead,imem)-Dat_All_0(ilon,ilat,it,16,ilead,imem) ! u200-u850
                     !else
                     !   X_inp(i,6)=Miss
                     !endif

                     X_inp(i,6)=Dat_All_0(ilon,ilat,it,16,ilead,imem)
                     ! 4 neighbor gridpoints
                     X_inp(i,7)=Dat_All_0(ilon-1,ilat,it,iprd,ilead,imem)          ! Model fcs
                     X_inp(i,8)=Dat_All_0(ilon+1,ilat,it,iprd,ilead,imem)          ! Model fcs
                     X_inp(i,9)=Dat_All_0(ilon,ilat-1,it,iprd,ilead,imem)          ! Model fcs
                     X_inp(i,10)=Dat_All_0(ilon,ilat+1,it,iprd,ilead,imem)          ! Model fcs

                     X_inp(i,11)=Dat_All_0(ilon-1,ilat,it,11,ilead,imem)   ! olr
                     X_inp(i,12)=Dat_All_0(ilon+1,ilat,it,11,ilead,imem)   ! olr
                     X_inp(i,13)=Dat_All_0(ilon,ilat-1,it,11,ilead,imem)   ! olr
                     X_inp(i,14)=Dat_All_0(ilon,ilat+1,it,11,ilead,imem)   ! olr

                     X_inp(i,15)=Dat_All_0(ilon-1,ilat,it,17,ilead,imem)   ! w700
                     X_inp(i,16)=Dat_All_0(ilon+1,ilat,it,17,ilead,imem)   ! w700
                     X_inp(i,17)=Dat_All_0(ilon,ilat-1,it,17,ilead,imem)   ! w700
                     X_inp(i,18)=Dat_All_0(ilon,ilat+1,it,17,ilead,imem)   ! 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"(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
         X_inp(1:Nsample1,j)=Inputs_X(j,1:Nsample1)
      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 Read_all_data (Dir,Model_nm)
      Character (Len=*) :: Dir,Model_nm
      Character (Len=255) :: Fname
      Real, Allocatable :: xlon(:),xlat(:),xtime(:)
      Integer :: imem,ilead,ivar,Nlo,Nla,Nt, k0, it1, it2

      k0=4    ! Bat dau chay mo hinh tu 4/2011
      Do imem=1,Nmem
         Do ilead=1,NLead
            write(Fname,"(A,I1.1,A,I2.2,A)") trim(Dir)//"lead",ilead,"_"//trim(Model_nm)//"_",imem,"_VN.nc"  ! clWRF_VN.nc"
            !it1=k0+ilead-1; it2=k0+7*Nmon+ilead
            Do ivar=1,Nvar
               Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,trim(vname(ivar)))
               if (imem==1.AND.ilead==1.AND.ivar==1) then
                  Ntime=Nyrs*Nmon ! Ntime=Nt
                  Allocate (Dat_All(Nlon,Nlat,Ntime,Nvar,NLead,Nmem)); Dat_All=Miss
               endif
               if (Nt/=Ntime) then; print*,"Something's wrong in Read_all_data"; stop; endif
               Dat_All(:,:,1:Nt,ivar,ilead,imem) = X3d(:,:,1:Nt)
               Deallocate (xlon,xlat,xtime,X3d)
            Enddo
         Enddo
      Enddo
   End Subroutine Read_all_data
   !--------------------------------------------------------------------------------
   Subroutine Data_standardization
      Integer :: ilon,ilat,ilead,it,ivar,imon, Nlead2, Nlo,Nla,nvr,nld
      Character (len = 15) :: lon_nm, lat_nm, lev_nm, t_nm
      Character (Len=2), Parameter :: mon(Nmon)=(/"01","02","03","04","05","06","07","08","09","10","11","12"/)
      Real, Allocatable :: xlon(:),xlat(:),xvar(:),xlead(:),A4d(:,:,:,:)
      Real    :: xdat,xmean,std,tmp

      Nlead2=Nlead+1     ! Lead=1-7 
      Allocate (Dat_mean(Nlon,Nlat,Nmon,Nvar,NLead2),Dat_std(Nlon,Nlat,Nmon,Nvar,NLead2))
      Fname=trim(Dir_Cli)//"Mean_and_Std_"//trim(Model_nm)//"_"//trim(yvarn)//".nc"
      Do imon=1,Nmon
         Call rd_ncdf_4d (Fname,xlon,xlat,xvar,xlead,Nlo,Nla,nvr,nld,A4d,"mean"//mon(imon))
         Dat_mean(:,:,imon,:,:)=A4d(:,:,:,:)
         Deallocate (xlon,xlat,xvar,xlead,A4d)
         Call rd_ncdf_4d (Fname,xlon,xlat,xvar,xlead,Nlo,Nla,nvr,nld,A4d,"std"//mon(imon))
         Dat_std(:,:,imon,:,:)=A4d(:,:,:,:)
         Deallocate (xlon,xlat,xvar,xlead,A4d)
      Enddo
      Allocate (Dat_All_0(Nlon,Nlat,Ntime,Nvar,NLead,Nmem)); Dat_All_0=Miss
      Do imem=1,Nmem
         Do ilead=1,NLead
            Do ivar=1,Nvar
               Do ilon=1,Nlon
                  Do ilat=1,Nlat
                     !xmean=Dat_mean(ilon,ilat,ivar,ilead)
                     !std  =Dat_std (ilon,ilat,ivar,ilead)
                     Do it=1,Ntime
                        imon=mod(it,Nmon); if (imon==0) imon=Nmon
                        xmean=Dat_mean(ilon,ilat,imon,ivar,ilead)
                        std  =Dat_std (ilon,ilat,imon,ivar,ilead)
                        tmp=Dat_All(ilon,ilat,it,ivar,ilead,imem)
                        if (tmp/=Miss.AND.xmean/=Miss.AND.std/=Miss) then
                           !xdat=(tmp - xmean)/std
                           xdat=(tmp-std)/(xmean-std)          ! Chuan hoa theo Max, Min
                           Dat_All_0(ilon,ilat,it,ivar,ilead,imem)=xdat
!if (ivar==1) print*,xmean,std,tmp,xdat
                        endif
                     Enddo
                  Enddo
               Enddo
            Enddo
!pause
         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"
   !-------------------------------------------------------------------------------------------