!
!  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   :: 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

   !--------------------------------------------------------------------------------
   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 
      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)") "/home/tanpv/Progs/562_VNU/test/NC/Out_"//trim(yvarn)//"_"//trim(Model_nm)//"_mon_",imon,"_lead_",ilead,".nc"
         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
            !print*,"ireg,imon,ilead=",ireg,imon,ilead
            Call Make_Data (trim(yvarn),ireg,imon,ilead)
            !write (save_file,"(3(A,I2.2),A)") "/home/tanpv/Progs/562_VNU/test/ANN_Save/Sav_"//trim(yvarn)//"_"//trim(Model_nm)//"_Reg_",ireg,"_mon_",imon,"_lead_",ilead,".txt"
            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)
            
            !mintmp=999999999.
            !Do i=2,7
            !   Nhidden=i
            !   Do k=600,3000,200
            !      Nepoc=k
            !      !Do k=600,3000,200
            !      Do j=1,1 !8
            !         learning_rate=0.001*real(j)
            !         learning_rate =   0.05
            !         Call ANN_2020
            !         if (Allocated(F1d)) Deallocate (F1d); if (Allocated(O1d)) Deallocate (O1d)
            !         Allocate (F1d(Nsample),O1d(Nsample))
            !         F1d(1:Nsample)=Yest(1:Nsample,1); O1d(1:Nsample)=Y_inp(1:Nsample) 
            !         tmp=RMAE(F1d,O1d,Nsample,Miss); 
            !         if (tmp<mintmp) then
            !            mintmp=tmp
            !            Op_Nhidden=Nhidden; Op_Nepoc=Nepoc; Op_learning_rate=learning_rate
            !         endif
            !         !print"(A,3I4)","ireg,imon,ilead=",ireg,imon,ilead
            !         print"(A,3I4,I3,F6.2,I6,F15.2)","ireg,imon,ilead,Nhidden,learning_rate,Nepoc,MAE=",ireg,imon,ilead,Nhidden,learning_rate,Nepoc,tmp
            !      Enddo
            !   Enddo
            !Enddo
            !print*, "Min RMAE for reg: ", ireg
            !print*,"Optimal values: "
            !print"(A,I3,F6.2,I8)","Nhidden,learning_rate,Nepoc=",Op_Nhidden,Op_learning_rate,Op_Nepoc

            Call Convert_to_lat_lon (ireg)
         Enddo
         Call Verification (ilead,imon)
         Call Write_Out (output_file,yvarn)
         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
   Fname=trim(Dir_Veri)//"Verification_scores_Training_"//trim(yvarn)//"_"//trim(Model_nm)//".nc" 
   Call Write_Verification (Fname)

   Open (3,file=trim(Dir_Veri)//"Score_Cor_AV_Training_"//trim(yvarn)//"_"//trim(Model_nm)//".txt",status="Unknown") 
   write (3,"(A)") "Percentage of grid points with positive correlation and added value"
   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),imon=1,Nmon)
         write (3,"(2A6,12F6.1)") "AV -  "," ",(Score_AV (ireg,imon,ilead),imon=1,Nmon)
      Enddo
   Enddo
   Close (3)
   !--------------------------------------------------------------------------------
      
   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
         !Do imem=1,Nz
            it=0
            Do iyr=1,Nyear
               Do imon=1,Nmon
                  it=it+1; Y4d(:,:,ivar,it)=X5d(:,:,imon,iyr,ivar)
               Enddo
            Enddo
         !Enddo
      Enddo
      !   var_longname=trim(var_ln)//v_ln(ivar)
      !   if (ivar==1) then
      !      Call wr_ncdf_3d (Fname,Lon,Lat,xtime,Nlon,Nlat,Ntime,Y3d,vname(ivar),t_unit, var_longname,"-", lon_nm, lat_nm, t_nm, Miss)
      Call wr_ncdf_4d (Fname,Lon,Lat,zlev,xtime,Nlon,Nlat,Nz,Ntime,Y4d,trim(yvarn), &
                      t_unit, var_ln, "-", lon_nm, lat_nm, z_nm, t_nm, Miss)
      !   else
      !      Call wr_ncdf_3d_apd (filenm,var,var_nm,nlon,nlat,nt, var_long_name, var_unit, lon_nm, lat_nm, t_nm, flag)
            !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 Write_Out (Fname,yvarn)
      Character (Len=*) :: Fname,yvarn
      Character (Len=150) :: vname, t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm
      Real, Allocatable, Dimension (:) :: zlev, xtime
      INteger :: i, Nz, Nyear

      Nz=3; Nyear=Nyrs !-1
      Allocate (zlev(Nz),xtime(Nyear))
      Do i=1,Nz; zlev(i)=i; Enddo; Do i=1,Nyear; xtime(i)=i+1981; Enddo
      t_unit="years since 1982-01-01 "; lon_nm="lon"; lat_nm="lat"; t_nm="year"; z_nm="var"
      var_ln="Outputs for "//trim(yvarn)//": 1) Observation; 2) Estimation; 3) Model forecast"
      Call wr_ncdf_4d (Fname,Lon,Lat,zlev,xtime,Nlon,Nlat,Nz,Nyear,X4d,yvarn, &
                       t_unit, var_ln, "-", lon_nm, lat_nm, z_nm, t_nm, Miss)
   End Subroutine Write_Out
   !--------------------------------------------------------------------------------
   Subroutine Verification (ilead,imon)
      Integer :: ilon,ilat,imon,ilead,ireg, ic0, ic1, ic2
      Integer, Parameter :: N=Nyrs !-1
      Real :: F1d(N),O1d(N),Model(N),Corr,Prob

      Do ilon=1,Nlon
         Do ilat=1,Nlat
            if (Int(Regs(ilon,ilat))>0) then
               F1d(1:N)=X4d(ilon,ilat,2,1:N); O1d(1:N)=X4d(ilon,ilat,1,1:N)  ! Observation
               Model(1:N)=X4d(ilon,ilat,3,1:N)       ! Model output
               Pr_ME  (ilon,ilat,ilead,imon) = ME  (F1d,O1d,N,Miss)
               Pr_MAE (ilon,ilat,ilead,imon) = MAE (F1d,O1d,N,Miss)
               Pr_RME (ilon,ilat,ilead,imon) = RME (F1d,O1d,N,Miss) 
               Pr_RMAE(ilon,ilat,ilead,imon) = RMAE(F1d,O1d,N,Miss) 
               Call Correlation (F1d,O1d,N,Corr,Prob,Miss)
               Pr_COR (ilon,ilat,ilead,imon) = Corr 
               Pr_Pval(ilon,ilat,ilead,imon) = Prob
               Pr_AV  (ilon,ilat,ilead,imon) = AV  (F1d,Model,O1d,N,Miss)
            endif
         Enddo
      Enddo 
      !
      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)>0) ic1=ic1+1
                  if (Pr_AV  (ilon,ilat,ilead,imon)>0) ic2=ic2+1
               endif
            Enddo
         Enddo
         Score_Cor(ireg,imon,ilead)=real(ic1)/real(ic0)*100.0
         Score_AV (ireg,imon,ilead)=real(ic2)/real(ic0)*100.0
      Enddo
   End Subroutine Verification
   !--------------------------------------------------------------------------------
   Subroutine Write_Verification (Fname)
      Character (Len=*) :: Fname
      Character (Len=150) :: vname, t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm
      Real, Allocatable, Dimension (:) :: zlev, xtime
      INteger :: i, Nz

      Allocate (zlev(NLead),xtime(Nmon))
      Do i=1,NLead; 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="Relative Mean Errors (RME)"
      Call wr_ncdf_4d (Fname,Lon,Lat,zlev,xtime,Nlon,Nlat,NLead,Nmon,Pr_RME,"rme", &
                       t_unit, var_ln, "-", lon_nm, lat_nm, z_nm, t_nm, Miss)
      var_ln="Relative Mean Absolute Errors (RMAE)"
      Call wr_ncdf_4d_apd (Fname,Pr_RMAE,"rmae", Nlon,Nlat,NLead,Nmon, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
      var_ln="Mean Errors (ME)"
      Call wr_ncdf_4d_apd (Fname,Pr_ME,"me", Nlon,Nlat,NLead,Nmon, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
      var_ln="Mean Absolute Errors (MAE)"
      Call wr_ncdf_4d_apd (Fname,Pr_MAE,"mae", Nlon,Nlat,NLead,Nmon, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
      var_ln="Correlation between corrected and observation"
      Call wr_ncdf_4d_apd (Fname,Pr_COR,"cor", Nlon,Nlat,NLead,Nmon, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
      var_ln="p-values of Correlation between corrected and observation"
      Call wr_ncdf_4d_apd (Fname,Pr_Pval,"pval", Nlon,Nlat,NLead,Nmon, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
      var_ln="Added values of bias correction in comparison to model output"
      Call wr_ncdf_4d_apd (Fname,Pr_AV,"av", Nlon,Nlat,NLead,Nmon, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
   End Subroutine Write_Verification
   !--------------------------------------------------------------------------------
   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!"
      !var_ln="Differences between model outputs and observarions used as climatological biases. z=7: Mean Obs."
      !Call wr_ncdf_4d (Fname,Lon,Lat,zlev,xtime,Nlon,Nlat,NLead+2,Nmon,Del_Cli,"bias", &
      !                 t_unit, var_ln, "-", lon_nm, lat_nm, z_nm, t_nm, Miss)
   End Subroutine Write_Bias_Cli
   !--------------------------------------------------------------------------------
   Subroutine Convert_to_lat_lon (ireg)
      Integer :: ireg,ilon,ilat,it,iyr

      it=0
      Do ilon=1,Nlon
         Do ilat=1,Nlat
            if (Int(Regs(ilon,ilat))==ireg) then
               Do iyr=1,Nyrs !-1
                  it=it+1
                  X4d(ilon,ilat,1,iyr)=Y_inp(it)
                  X4d(ilon,ilat,2,iyr)=Yest(it,1)
                  X4d(ilon,ilat,3,iyr)=X_inp(it,1)
               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
print*,"Beginning"
      Do it=1,3             !iyr=1,Nyear
         Do iyr=1,Nyear     ! imem=1,Nmem
            !Do imem=1,Nmem  !it=1,3
               X5d(:,:,imon,iyr,it)=X4d(:,:,it,iyr)
            !Enddo
            !X6d(:,:,imon,iyr,Nmem+1,it)=X4d(:,:,it,iyr)
         Enddo
      Enddo
      ! Climatological Bias Correction:
      Do ilon=1,Nlon
         Do ilat=1,Nlat
            if (Int(Regs(ilon,ilat))>0) then
               if (X4d(ilon,ilat,3,iyr)/=Miss.AND.Bias_Cli(ilon,ilat,ilead,imon)/=Miss) then
                  X5d(ilon,ilat,imon,iyr,4)=X4d(ilon,ilat,3,iyr)-Bias_Cli(ilon,ilat,ilead,imon) ! X1d = Model-Bias = Corrected
                  if (X5d(ilon,ilat,imon,iyr,4)<0) X5d(ilon,ilat,imon,iyr,4)=0     ! rainfall must be not negative
               else
                  X5d(ilon,ilat,imon,iyr,4)=Miss 
               endif
            endif
         Enddo
      Enddo
print*,"End"
   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=8                 ! 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(ilon,ilat,it,ipredictand,ilead)   ! Observation
                     X_inp(i,1)=Dat_All(ilon,ilat,it,iprd,ilead)          ! Model fcs
                     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
      !i=0
      !Do it=1,Nspl
      !   if (Y_inp(it)/=Miss.AND..not.Miss_Ar1(X_inp(it,1:Mx),Mx,Miss)) i=i+1
      !Enddo
      !Nsample1=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
         !if (Y_inp(it)/=Miss.AND..not.Miss_Ar1(X_inp(it,:),Mx,Miss)) then
            i=i+1
            Outputs_Y(1,i)=Y_inp(it)
            Do j=1,Mx
               Inputs_X(j,i)=X_inp(it,j)
            Enddo
         !endif
      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

      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:Nt0)=Dat_All(ilon,ilat,1:Nt0,iob,ilead)
                  F(1:Nt0)=Dat_All(ilon,ilat,1:Nt0,ifc,ilead)
                  Do imon=1,Nmon
                     iyr=0; 
                     Do it=1,Nt0
                        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)
                     Del_Cli(ilon,ilat,ilead,imon)=MeanF-MeanO
                     if (ilead==1) Del_Cli(ilon,ilat,NLead+2,imon)=MeanO
                  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

      !Dir="/work/users/hmo/tanpv/From_HOME/562/Model_Data/"
      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(:,:,:,ivar,ilead) = X3d(:,:,:)
            Dat_All(:,:,1:Nt,ivar,ilead) = X3d(:,:,1:Nt)
            Deallocate (xlon,xlat,xtime,X3d)
         Enddo
      Enddo
   End Subroutine Read_all_data
   !--------------------------------------------------------------------------------
   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/Stdev.f90"
   Include "/home/tanpv/Mylibs/Ave.f90"
   !--------------------------------------------------------------------------------
   END
   Include "/home/tanpv/Mylibs/Student_Prob.f90"
   !-------------------------------------------------------------------------------------------