!
!  CHUONG TRINH Du bao So luong XTND bang cac Mo hinh Thong ke
!  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: La cac EOF (3 EOF dau tien) cua 14 ClimIDX. So nhan to du
!  tuyen = 3 thang truoc thang du bao x 3 EOF = 9 nhan to
!  Cac nhan to du bao thuc te se duoc chon tu 9 nhan to du tuyen boi HQTB

   PROGRAM TC_Fcs_2020
   Use Netcdf

   Integer, Parameter :: Nmon=12, Nfile=Nmon, Nmodels=9, MMax=10
   Integer            :: Yr1, Yr2, Fcs_mon, Nyr1, NMax, Ncase, Ntc_var
   Real, Parameter    :: Miss=-99.

   Integer, Allocatable :: Yr01(:), Yr02(:)

   REAL, ALLOCATABLE :: X(:,:), X1(:)              ! Input Data Arrays
   INTEGER, ALLOCATABLE :: IDX(:)                  ! Index of variables for stepwise regression

   Character (Len=70) :: FINP !, FOUT                ! Filename for input and output

   REAL, ALLOCATABLE :: Xorg(:,:), Yorg(:), Yhq(:)
   Character (Len=8), Allocatable :: Ten(:), Varname(:)
   REAL, ALLOCATABLE :: Beta_tmp(:)

   Character (Len=70) :: TIEUDE
   Character (Len=2)  :: ST1,ST2
   Character (Len=20) :: FMT

   !--------------------------------------------------------------------------------

   Integer :: ifile                    ! So file thuc te
   Integer :: Mvars, N, Mbien              ! So bien doc lap, Dung luong mau va Tat ca cac Bien (Mvars+1)

   Real, Allocatable :: Beta(:,:,:)   ! Cac he so HQ cua 3 mo hinh HQTT, LAD, MINIMAX,
   Real, Allocatable :: Yout(:,:,:)   ! Outputs of Yorg and Y-estimated by 4 models
   Real, Allocatable, Dimension(:,:) :: ME, MAE, RMSE, RME, Corre, Pvalue  ! Verification

   Integer :: Nyear, imon      ! Tong so nam co SL va chi so thang
   !--------------------------------------------------------------------------------
   ! ANN:
   character (len=250) :: input_file, output_file, save_file
   logical :: training_mode

   ! Use for calculating scaling factors
   Integer :: Nsample, Nepoc

   Real,allocatable, Dimension(:,:) :: Inputs, Outputs
   Real,allocatable, Dimension(:)   :: min_in, max_in, min_out, max_out
   Real :: mean_error
   !-------------------------------------------------------------------------------
    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_rate,momentum, error
    logical :: debug
    Character (Len=250) :: Y_out_Fname_NC, Y_out_Fname_TXT, Beta_Fname, Veri_Fname_NC, Veri_Fname_TXT
 
   !--------------------------------------------------------------------------------

   Real, Allocatable :: Mean_Obs(:), Std_Obs(:)
   
   Character (Len=5), Parameter :: Mdel_nm(Nmodels+1)=(/ &
             "Obs","MLR","LAD","LMV","ANN","E12","E123","E124","E14","E1234"/)

   Integer :: Yr_Sta(Nmon), Nyrs(Nmon), NVars(Nmon)   ! Nam bat dau, Dung luong mau, So bien DL 
   Character (Len=250) :: Dir, Dir_fc, Namelist_file
   Character (Len=5)   :: tc_var
   Character (Len=7)   :: x_var
   Character (Len=5), Allocatable   :: tc_varn(:)  !="ntc3m"
   Character (Len=7), Allocatable   :: x_varn(:)   !="eof"
   Character (Len=250), Allocatable :: Dir_Inp(:), Dir_Fcs(:)

   Character (Len=250) :: Fnames(Nfile), Fnm
   Character (Len=250) :: FnameOut(Nfile), FnameSave(Nfile), FnameInp(Nfile)

   Namelist /Const/ Ncase, Ntc_var
   Namelist /Arrays/ Yr01, Yr02, Fcs_mon, tc_varn, x_varn, Dir_Inp, Dir_Fcs

   Integer :: itc, jtc
   !--------------------------------------------------------------------------------
   ! For ensemble:
   Integer, Parameter :: Model=4 !Nmodels !4            ! MLR, LAD, LMV, ANN
   Integer :: Yrm, Yrx, Nyrx, iy1,iy2,iy, Nens, Nens1, Del   !
   Real, Allocatable :: TC_Obs(:,:,:,:)        ! (Year, Mon, Case, TC_Cha)
   Real, Allocatable :: TC_Fcs(:,:,:,:,:)   ! (Year, Model, Case, Mon, TC_Cha) - Model=4
   Real, Allocatable :: TC_Obs_Fcs(:,:,:,:) ! (Year, Model*Case, Mon, TC_Cha) 
   Real, Allocatable :: X4d(:,:,:,:), Xtc(:), Xmem(:), Xyr(:), Xmon(:), Xtc_var(:)
   !--------------------------------------------------------------------------------

   ! Read Namelist:
   Namelist_file="Namelist_Forecast_Ope.txt" !Namelist_Forecast.txt"
   Open (1,file=Namelist_file,Status="Old")

   Read (1, Const)     ! Get Ncase, Ntc_var
   Allocate (Yr01(Ncase),Yr02(Ncase),tc_varn(Ntc_var),x_varn(Ncase),Dir_Inp(Ncase),Dir_Fcs(Ncase))
   Read (1, Arrays)

   !--------------------------------------------------------------------------------
   ! For Probability Fcs:
   Yrm=9999; Yrx=0
   Do icase=1,Ncase
      if (Yr01(icase)<Yrm) Yrm=Yr01(icase)
      if (Yr02(icase)>Yrx) Yrx=Yr02(icase)
   Enddo
   Nyrx = Yrx-Yrm+1
   Nens = Model*Ncase+1   ! Including Obs
   Allocate (TC_Obs(Nyrx,Nmon,Ncase,Ntc_var), TC_Fcs(Nyrx,Nmon,Ncase,Model,Ntc_var))
   Allocate (TC_Obs_Fcs(Nyrx,Nmon,Nens,Ntc_var))
   TC_Obs = Miss; TC_Fcs = Miss; TC_Obs_Fcs = Miss
   !------------------- For ANN -------------------------------------------------------------
     Nhidden       =     2        ! Number of hidden layers
     learning_rate =   0.5
     momentum      =     0
     training_mode =.false.       ! For training
     Nepoc         =  8000        ! Iterative training times
   !--------------------------------------------------------------------------------

    Do icase=1,Ncase
       Yr1=Yr01(icase); Yr2=Yr02(icase)
       NMax=Yr2-Yr1+1
       Allocate (Xorg(NMAX,MMax), Yorg(NMAX), Yhq(NMAX))
       Allocate (Ten(MMAX), Varname(MMax))
       Allocate (Beta_tmp(MMax+1))

       Dir=trim(Dir_Inp(icase)); x_var=trim(x_varn(icase))
       Dir_fc=trim(Dir_Fcs(icase))

       Do jtc=1,Ntc_var
          tc_var=trim(tc_varn(jtc))

          Do i=1,Nfile
             write (Fnames(i),"(A,A,I2.2,A)") trim(tc_var),"_",i,".txt"
          Enddo
          Do i=1,Nfile 
             FnameInp(i)=trim(Dir_fc)//trim(Fnames(i)); FnameOut(i)=trim(Dir)//"Out_"//trim(Fnames(i)); FnameSave(i)=trim(Dir)//"Sav_"//trim(Fnames(i))
          Enddo
          Y_out_Fname_NC   = "Dat_Ope/Yout_fc_"//trim(tc_var)//"_"//trim(x_var)//".nc"
          Y_out_Fname_TXT  = "Dat_Ope/TC_out_fc_"//trim(tc_var)//"_"//trim(x_var)//".txt"
          Beta_Fname       = "Dat_Ope/Beta_"//trim(tc_var)//"_"//trim(x_var)//".nc"
          Veri_Fname_NC    = "Dat_Ope/Veri_Results_fc_"//trim(tc_var)//"_"//trim(x_var)//".nc"
          Veri_Fname_TXT   = "Dat_Ope/Veri_Results_fc_"//trim(tc_var)//"_"//trim(x_var)//".txt"
          !--------------------------------------------------------------------------------

          Allocate (Mean_Obs(Nfile), Std_Obs(Nfile))
          Allocate (Yout(Nfile,Nmax,0:Nmodels))
          Allocate (ME(Nfile,Nmodels),MAE(Nfile,Nmodels),RMSE(Nfile,Nmodels))
          Allocate (RME(Nfile,Nmodels),Corre(Nfile,Nmodels),Pvalue(Nfile,Nmodels))     

          Nyear=0; Yout(:,:,:) = Miss;
          Call Read_Beta (trim(Beta_Fname))

          Do ifile=Fcs_mon,Fcs_mon !1, Nfile

             Call Read_Data (ifile, FnameInp(ifile)) 
             imon = ifile; Mvars = NVars(imon)-1; Mbien=NVars(imon); N=Nyrs(imon)

             Yout(ifile,1:N,0) = Yorg(1:N) 
             Mean_Obs(ifile) = Ave  (Yorg(1:N), N, Miss)
             Std_Obs(ifile)  = Stdev(Yorg(1:N), N, Miss) 
 
             print*,"Calling Calc_HQTT..."
             Call Calc_HQTT (ifile)        ! HQTT
             Yout(ifile,1:N,1) = Yhq(1:N)
             print*,"Calling Calc_LAD..."

             Call Calc_LAD (ifile)        ! Least Absolute Deviation
             Yout(ifile,1:N,2) = Yhq(1:N)
             print*,"Calling Calc_LMV..."

             Call Calc_LMV (ifile)        ! Minimax
             Yout(ifile,1:N,3) = Yhq(1:N)
             !--------------------------ANN-------------------------
             debug = .false.
             input_file  = trim(FnameInp(ifile))
             output_file = trim(FnameOut(ifile))
             save_file   = trim(FnameSave(ifile))

             print*, "Calculating ",input_file
             Call controller_init

             if (training_mode) then       ! .T. = Training, .F. = forecast (Define from Namelist)
                Call controller_train_data
             endif

             write(*,*) "Writting output..."
             Call controller_output_data !(Yhq)
             Call controller_deinit
             Yout(ifile,1:N,4) = Yhq(1:N)
             !--------------------------ANN-------------------------

          Enddo

          ! Set all estimated negative TC-counts to Zero:
          Where (Yout/=Miss.AND.Yout<0) Yout=0

           ! Small ensemble:
          Do ifile=Fcs_mon,Fcs_mon  !1, Nfile
             Nyear=Nyrs(ifile)
             Do iyr=1, Nyear
                if (Yout(ifile,iyr,1)/=Miss.AND.Yout(ifile,iyr,2)/=Miss) &
                   Yout(ifile,iyr,5)=(Yout(ifile,iyr,1)+Yout(ifile,iyr,2))/2.
 
                if (Yout(ifile,iyr,1)/=Miss.AND.Yout(ifile,iyr,2)/=Miss.AND.Yout(ifile,iyr,3)/=Miss) &
                   Yout(ifile,iyr,6)=(Yout(ifile,iyr,1)+Yout(ifile,iyr,2)+Yout(ifile,iyr,3))/3.

                if (Yout(ifile,iyr,1)/=Miss.AND.Yout(ifile,iyr,2)/=Miss.AND.Yout(ifile,iyr,4)/=Miss) &
                   Yout(ifile,iyr,7)=(Yout(ifile,iyr,1)+Yout(ifile,iyr,2)+Yout(ifile,iyr,4))/3.

                if (Yout(ifile,iyr,1)/=Miss.AND.Yout(ifile,iyr,4)/=Miss) &
                   Yout(ifile,iyr,8)=(Yout(ifile,iyr,1)+Yout(ifile,iyr,4))/2.
        
                if (Yout(ifile,iyr,1)/=Miss.AND.Yout(ifile,iyr,2)/=Miss.AND.Yout(ifile,iyr,3)/=Miss.AND.Yout(ifile,iyr,4)/=Miss) &
                   Yout(ifile,iyr,9)=(Yout(ifile,iyr,1)+Yout(ifile,iyr,2)+Yout(ifile,iyr,3)+Yout(ifile,iyr,4))/4.
             Enddo

             ! Tinh sai so va HSTQ giua Yorg va Yhq
             Yorg(1:Nyear) = Yout(ifile,1:Nyear,0) 
             Do imod=1, Nmodels  
                print*,"Calling Verification of Model: ", imod
                Yhq(1:Nyear) = Yout(ifile,1:Nyear,imod) ! ; Where (Yhq/=Miss.AND.Yhq<0) Yhq=0
             Enddo
          Enddo

           ! For all cases including ensemble
           Do ifile=Fcs_mon,Fcs_mon  ! 1, Nfile
              Nyear=Nyrs(ifile)
              Do iyr=1,Nyear !Yrm,Yrx
                 iy1=Yr1+iyr-1
                 Do iy=1,Nyrx
                    iy2=Yrm+iy-1
                    !print*,"iyr,iy,iy1,iy2=",iyr,iy,iy1,iy2
                    if (iy2==iy1) then
                       TC_Obs(iy,ifile,icase,jtc)        =Yout(ifile,iyr,0)
                       TC_Fcs(iy,ifile,icase,1:Model,jtc)=Yout(ifile,iyr,1:Model)
                       !print*,"iyr,iy,iy1,iy2=",iyr,iy,iy1,iy2
                       exit
                    endif
                 Enddo
              Enddo
           Enddo
 

          print*,"Writing to TXT file ..."
          Call Write_Yorg_Yhq_to_TXT (trim(Y_out_Fname_TXT)) !"TC_out.txt")

          Deallocate (Mean_Obs,Std_Obs)
          Deallocate (Beta,Yout,ME,MAE,RMSE,RME,Corre,Pvalue)

       Enddo

       Deallocate (Xorg,Yorg,Yhq,Ten,Varname,Beta_tmp)

    Enddo

     ! Write out for probability fcs:
     ! TC_Obs_Fcs(Nyrx,Nmon,Nens,Ntc_var)

     Nens1 = Nens+1
     Allocate (X4d(Nens1,Nmon,Ntc_var,Nyrx), Xmem(Nens1), Xyr(Nyrx), Xmon(Nmon), Xtc_var(Ntc_var),Xtc(Nens-1))

     Open (3,file="Dat_Ope/TC_nTC_Obs_Fcs_Forecast_All.txt",status="Unknown")
     write (3,*) "TC/Ntc obs and fcs"
     Do jtc=1,Ntc_var
        if (mod(jtc,2)==0) then; k_mon=5; else; k_mon=2; endif
        if (jtc<=2) then; Del=1; else; Del=4; endif       ! Khoang Del tinh XS
        Do imon=Fcs_mon,Fcs_mon  ! 1,Nmon
           write(3,"(3A,I2.2,A,I2.2,A,I2.2)") "TCvarname: ",tc_varn(jtc)," Month: ",imon,"->",imon,"+",k_mon
           write(3,"(A4,A6,2x,36(A4,I2.2,1x))") "Year","Obs",(("Mem_",iy1),iy1=1,Nens-1)
           Do iyr=1,Nyrx
              Do icase=1,Ncase
                 if (TC_Obs(iyr,imon,icase,jtc)>=0) then
                    TC_Obs_Fcs(iyr,imon,1,jtc) = TC_Obs(iyr,imon,icase,jtc); exit
                 endif
              Enddo
              iy=1
              Do icase=1,Ncase
                 Do imod=1,Model
                    iy=iy+1
                    TC_Obs_Fcs(iyr,imon,iy,jtc) = TC_Fcs(iyr,imon,icase,imod,jtc)
                 Enddo
              Enddo
              X4d(1:Nens,imon,jtc,iyr)=TC_Obs_Fcs(iyr,imon,1:Nens,jtc)
              Xtc(1:Nens-1)=TC_Obs_Fcs(iyr,imon,2:Nens,jtc)
              X4d(Nens+1,imon,jtc,iyr)=TC_prob (Xtc,Nens-1,Del,Miss)
              write(3,"(I4,F6.1,2x,36(F6.1,1x))") iyr+Yrm-1,TC_Obs_Fcs(iyr,imon,1:Nens,jtc)

              !print"(I4,5F6.1,2x,20F6.1)",iyr+Yrm-1,TC_Obs(iyr,imon,1:Ncase,jtc),TC_Fcs(iyr,imon,1:Ncase,1:Model,jtc)
              !print"(I4,F6.1,2x,20F6.1)", iyr+Yrm-1,TC_Obs_Fcs(iyr,imon,1:Nens,jtc)
           Enddo
           !print*,"Mon:",imon," tcvar: ",tc_varn(jtc); pause
        Enddo
     Enddo
     Close (3)

     Fnm="Dat_Ope/TC_nTC_Obs_Fcs_Forecast_All.nc"
     Do jtc=1,Ntc_var; Xtc_var(jtc)=jtc; Enddo
     Do imon=1,Nmon; Xmon(imon)=imon; Enddo
     Do iyr=1,Nyrx; Xyr(iyr)=iyr+Yrm-1; Enddo
     Do iy=1,Nens+1; Xmem(iy)=iy; Enddo
     Call wr_ncdf_4d (Fnm,Xmem,Xmon,Xtc_var,Xyr,Nens+1,Nmon,Ntc_var,Nyrx,X4d,"tc", &
                       "years since 2010-01-01 00","TC/NTC obs and forecast","-","mem","mon","tcvar","yr",Miss)
   !--------------------------------------------------------------------------------
      
   CONTAINS
   !--------------------------------------------------------------------------------
   Function TC_prob (X,N,Del,Miss)
      Integer :: N, i,j, ic, Xmax, Xmin, Del, Nbin, tmp
      Real :: X(N), Miss, Xtc
      Integer, Allocatable :: Bin(:)
      Real, Allocatable    :: Xbin(:)

      Xmax=-9999; Xmin=9999
      Do i=1,N
         if (X(i)==Miss) Cycle
         tmp=Nint(X(i))
         if (tmp>Xmax) Xmax=tmp; if (tmp<Xmin) Xmin=tmp
      Enddo
      Nbin=(Xmax-Xmin)/Del+2
      
      Allocate (Bin(Nbin),Xbin(Nbin-1))
      Do j=1,Nbin; Bin(j)=Xmin+(j-1)*Del; Enddo
      Xbin(:)=0
      Do j=1,Nbin-1
         ic=0
         Do i=1,N
            if (X(i)/=Miss) then
               ic=ic+1; tmp=Nint(X(i))
               if (Del==1) then
                  if (tmp==Bin(j)) Xbin(j)=Xbin(j)+1
               else
                  if (tmp>=Bin(j).AND.tmp<Bin(j+1)) Xbin(j)=Xbin(j)+1
               endif
            endif
         Enddo
         Xbin(j)=Xbin(j)/real(ic)
      Enddo
      Xtc=-9999
      Do i=1,Nbin
         if (Xbin(i)>Xtc) then; Xtc=Xbin(i); j=i; endif
      Enddo
      if (Del==1) then
         TC_prob=Bin(j)
      else
         if (j<Nbin) then 
            TC_prob=(Bin(j)+Bin(j+1))/2.
         else
            TC_prob=(Bin(j))
         endif
      endif
   End Function TC_prob
   !--------------------------------------------------------------------------------
   Subroutine Read_Beta (Fname)
     Character (Len=*) :: Fname
     Real, Allocatable :: Xlon(:), Xlat(:), Xtime(:)
     Integer :: Nfi, Ncoef

     Call rd_ncdf_3d (Fname,Xlon,Xlat,Xtime, Nfi, N_model, Ncoef, Beta, "Beta")
     print*, Nfi, N_model, Ncoef, Beta
     Deallocate (Xlon,Xlat,Xtime)
   End Subroutine Read_Beta
   !--------------------------------------------------------------------------------

Subroutine Write_Yorg_Yhq_to_TXT (Fname)
  Character (Len=*) :: Fname
  Real, Allocatable :: X(:,:,:) ! (Nfile, Nmodels+1, N)
  Integer :: ifile, iyr, imon, imod, it1, it2, Nt
  ! Yr_Sta(Nmon), Nyrs(Nmon)
  
  Nyear=Nmax
  Open (3,file=trim(Fname),status="Unknown")
  Allocate (X(1:Nfile,1:Nmodels+1,1:Nyear)); X(:,:,:) = Miss

  Do ifile=Fcs_mon,Fcs_mon  ! 1, Nfile

     write(3,*) "Input file: ",trim(FnameInp(ifile))
     write(3,"(A5,10A7)") "Yrs", (Mdel_nm(imod),imod=1,Nmodels+1)
     imon=ifile
     it1=Yr_Sta(imon)-Yr1+1; it2=it1+Nyrs(imon)-1; Nt=Nyrs(imon)
     Do imod=1, Nmodels+1
        !X(ifile,imod,1:Nyear) = Yout(ifile,1:Nyear,imod-1)
        X(ifile,imod,it1:it2) = Yout(ifile,1:Nt,imod-1)
     Enddo

     Do iyr=1, Nyear
        write(3,"(I5,10F7.1)") iyr+Yr1-1, (X(ifile,imod,iyr), imod=1,Nmodels+1)
     Enddo
  Enddo
  Deallocate (X); Close (3)
  print*,"Write_Yorg_Yhq_to_TXT... Done"
End Subroutine Write_Yorg_Yhq_to_TXT
!--------------------------------------------------------------------------------
Subroutine Write_Yorg_Yhq_to_NC (Fname)
  Character (Len=255) :: t_unit, var_ln
  Character (Len=*) :: Fname
  Real, Allocatable :: X(:,:,:) ! (Nfile, Nmodels+1, N)
  Real, Allocatable :: Xlon(:), Xlat(:), Xtime(:)
  Integer :: ilon, ilat, it, imon, it1, it2, Nt

  ! Yr_Sta(Nmon), Nyrs(Nmon)
  Nyear=Nmax
  Allocate (X(1:Nfile,1:Nmodels+1,1:Nyear),Xlon(Nfile),Xlat(Nmodels+1),Xtime(Nyear)); X(:,:,:) = Miss
  Do ilon=1, Nfile; Xlon(ilon)=ilon; Enddo     ! So truong hop du bao: 12 thang x 2 vung = 24
  Do ilat=1, Nmodels+1; Xlat(ilat)=ilat; Enddo ! Yorg + 4 models
  Do it=1, Nyear; Xtime(it)=Yr1+it-1; Enddo              ! So nam du bao = 1962..2010
  write (t_unit,"(A,I4.4,A)") "years since ",Yr1,"-01-01 "
  var_ln="Observed TC-count and its Estimation by different models (MLR, LAD, LMV, ANN)"
  X=Miss
  Do ifile=1, Nfile
     imon=ifile !mod(ifile,12); if (imon==0) imon=12 
     it1=Yr_Sta(imon)-Yr1+1; it2=it1+Nyrs(imon)-1; Nt=Nyrs(imon) 
     print*,"Write_Yorg_Yhq_to_NC: it1,it2,Nt,Nyear:",it1,it2,Nt,Nyear
     Do imod=1, Nmodels+1
        !X(ifile,imod,1:Nyear) = Yout(ifile,1:Nyear,imod-1)
        X(ifile,imod,it1:it2) = Yout(ifile,1:Nt,imod-1)
     Enddo
  Enddo
  Call wr_ncdf_3d(trim(Fname),XLon,XLat,Xtime, Nfile,Nmodels+1,Nyear, X, "tc", trim(t_unit), trim(var_ln),"1","mon","model","year",Miss) 
  ! wr_ncdf_3d (filenm,lon,lat,t,nlon,nlat,nt,var,var_nm,t_unit, var_long_name, var_unit, lon_nm, lat_nm, t_nm, flag)
  Deallocate (Xlon, Xlat, Xtime, X)
End Subroutine Write_Yorg_Yhq_to_NC 
!--------------------------------------------------------------------------------

Subroutine Calc_HQTT (ifile)
  Beta_tmp(1:Mvars+1) = Beta(ifile,1,1:Mvars+1)
  print*,Beta_tmp(1:Mvars+1)
  Do i=1, N
     Yhq(i) = 0
     Do j=2, Mvars+1
        Yhq(i) = Yhq(i)+Beta_tmp(j)*Xorg(i,j-1) ! Beta_tmp(2..Mvars+1)=A1,...Am
     Enddo
     Yhq(i) = Yhq(i)+Beta_tmp(1)                ! Beta_tmp(1)=A0
  Enddo
End Subroutine Calc_HQTT
!--------------------------------------------------------------------------------
Subroutine Calc_LAD (ifile)
  Integer :: i, j, IOPT

  Beta_tmp(1:Mvars+1) = Beta(ifile,2,1:Mvars+1)
  Do i=1, N
     Yhq(i) = 0
     Do j=2, Mvars+1
        Yhq(i) = Yhq(i)+Beta_tmp(j)*Xorg(i,j-1) ! Beta_tmp(2..Mvars+1)=A1,...Am
     Enddo
     Yhq(i) = Yhq(i)+Beta_tmp(1)                ! Beta_tmp(1)=A0
  Enddo
End Subroutine Calc_LAD
!--------------------------------------------------------------------------------
Subroutine Calc_LMV (ifile)
  Integer :: i, j, IOPT

  Beta_tmp(1:Mvars+1) = Beta(ifile,3,1:Mvars+1)
  Do i=1, N
     Yhq(i) = 0
     Do j=2, Mvars+1
        Yhq(i) = Yhq(i)+Beta_tmp(j)*Xorg(i,j-1) ! Beta_tmp(2..Mvars+1)=A1,...Am
     Enddo
     Yhq(i) = Yhq(i)+Beta_tmp(1)                ! Beta_tmp(1)=A0
  Enddo
End Subroutine Calc_LMV
!--------------------------------------------------------------------------------

Subroutine Read_Data (ifile, Finp) ! Calc_Stepwise_Regression (ifile, Finp, Fout)
   Character (Len=*) :: Finp !, Fout
   Integer :: ifile, iv_y, Nc, Nyr1   ! Chi so bien phu thuoc
   Real, Allocatable :: Xtmp(:)
   Real              :: Tmp

      OPEN(1,FILE=trim(FINP),Status="Old")

      PRINT*," READ IN DATA FROM FILE "//trim(FINP) !,"  Output file:",trim(FOUT)

      READ (1,"(A)") TIEUDE
      Read(1,*) Nyr1, N, Mbien
      ALLOCATE (IDX(Mbien), Xtmp(Mbien-1))

      Read(1,*) (Idx(i),i=1,Mbien)
      Read (1,"(20A10)") (Ten(i),i=1,Mbien)
      Nc=0
   L1:   DO I=1,N
         READ(1,*) Tmp, (Xtmp(j), j=1,Mbien-1)
         !if (Tmp==Miss) Cycle L1
         Do j=1, Mbien-1
            if (Xtmp(j)==Miss) Cycle L1
         Enddo
         Nc=Nc+1;
         Yorg(Nc)=Tmp; Xorg(Nc,1:Mbien-1)=Xtmp(1:Mbien-1)
!         READ(1,*) Yorg(i), (Xorg(i,j), j=1,Mbien-1) ! (X(I,J),J=1,Mbien)
      ENDDO L1
      N = Nc
      CLOSE(1)
      Deallocate (IDX,Xtmp)
      !
      Yr_Sta(ifile)=Nyr1; Nyrs(ifile)=N; NVars(ifile)=Mbien

End Subroutine Read_Data 
!--------------------------------------------------------------------------------
     
Include "SUB/Sub_HQTT.f90"
Include "SUB/Sub_LAD.f90"
Include "SUB/Sub_LMV.f90"
Include "SUB/Calc_Correlation_Coefficient.f90"
Include "SUB/sub_wr_ncdf_Tan.F90" 
Include "SUB/sub_rd_ncdf_Tan.F90"
Include "SUB/Sub_neural.f90"
Include "SUB/Stdev.f90"
Include "SUB/Ave.f90"
!--------------------------------------------------------------------------------
END
Include "SUB/Student_Prob.f90"
!-------------------------------------------------------------------------------------------