!
!  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, 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=200) :: 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 

   Integer :: icase 
   Character (Len=250) :: Dir, 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(:)

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

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

   Integer :: itc, jtc, k_mon
   !--------------------------------------------------------------------------------
   ! For ensemble:
   Integer, Parameter :: Model=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_Training_Ope.txt" !Namelist_Training.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))
   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
   !print*,"Yrm,Yrx,Nyrx:",Yrm,Yrx,Nyrx; pause
   !------------------- For ANN -------------------------------------------------------------
     Nhidden       =     2        ! Number of hidden layers
     learning_rate =   0.5
     momentum      =     0
     training_mode = .true.       ! 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))

        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)//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_"//trim(tc_var)//"_"//trim(x_var)//".nc"
           Y_out_Fname_TXT  = "Dat_Ope/TC_out_"//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_"//trim(tc_var)//"_"//trim(x_var)//".nc"
           Veri_Fname_TXT   = "Dat_Ope/Veri_Results_"//trim(tc_var)//"_"//trim(x_var)//".txt"
           !--------------------------------------------------------------------------------

           Allocate (Mean_Obs(Nfile), Std_Obs(Nfile))
           Allocate (Beta(Nfile,Nmodels,0:Mmax), 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;

           Do ifile=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        ! HQTT
              Beta(ifile,1,0) = Beta_tmp(1); Beta(ifile,1,1:Mvars) = Beta_tmp(2:Mvars+1)
              Yout(ifile,1:N,1) = Yhq(1:N)
              print*,"Calling Calc_LAD..."

              Call Calc_LAD         ! Least Absolute Deviation
              Beta(ifile,2,0) = Beta_tmp(Mvars+1); Beta(ifile,2,1:Mvars) = Beta_tmp(1:Mvars)
              Yout(ifile,1:N,2) = Yhq(1:N)
              print*,"Calling Calc_LMV..."

              Call Calc_LMV         ! Minimax
              Beta(ifile,3,0) = Beta_tmp(Mvars+1); Beta(ifile,3,1:Mvars) = Beta_tmp(1:Mvars)
              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

           ! For probability Fcs:
           ! TC_Obs(Nyrx,Nmon,Ncase,Ntc_var), TC_Fcs(Nyrx,Nmon,Ncase,Model,Ntc_var) 
           Do ifile=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
           
           ! Small ensemble:
           Do ifile=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
                 Call Calc_Verification (Yorg, Yhq, Nyear, ME(ifile,imod), MAE(ifile,imod), RMSE(ifile,imod), RME(ifile,imod), Corre(ifile,imod), Pvalue(ifile,imod))
              Enddo
           Enddo

           print*,"Write_Yorg_Yhq_to_NC..."
           Call Write_Yorg_Yhq_to_NC (trim(Y_out_Fname_NC)) ! "Yout.nc")                      ! Gia tri TC observed va estimated

           print*,"Write Beta..."
           Call Write_Model_Coeffs_to_NC (trim(Beta_Fname)) !"Beta.nc")

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

           print*,"Write_Verification_Results_to_NC"
           Call Write_Verification_Results_to_NC (trim(Veri_Fname_NC)) !"Veri_Results.nc")  ! Cac chi so danh gia
  
           Call Write_Verification_to_TXT (trim(Veri_Fname_TXT)) !"Veri_Results.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(Nens,Nmon,Ntc_var,Nyrx), Xtc(Nyrx), Xmem(Nens), Xyr(Nyrx), Xmon(Nmon), Xtc_var(Ntc_var))
     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_Training_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=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,20(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
              write(3,"(I4,F6.1,2x,20(F6.1,1x))") iyr+Yrm-1,TC_Obs_Fcs(iyr,imon,1:Nens,jtc)
              !X4d(1:Nens,imon,jtc,iyr)=TC_Obs_Fcs(iyr,imon,1:Nens,jtc)
              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)

              !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_Training_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; Xmem(iy)=iy; 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 1980-01-01 00","TC/NTC obs and forecast","-","mem","mon","tcvar","yr",Miss)
     !Allocate (X4d(Nens,Nmon,Ntc_var,Nyrx), Xtc(Nyrx), Xmem(Nens), Xyr(Nyrx), Xmon(Nmon), Xtc_var(Ntc_var))
   !--------------------------------------------------------------------------------
      
   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 Write_Verification_to_TXT (Fname)
  Integer, Parameter :: Nmon=12
  Character (Len=*) :: Fname
  Character (Len=50) :: Region ! (/"VES_N15","VES_S15"/) 
  Integer :: imod, ifile, imon
  
  Open (3,file=trim(Fname), Status="Unknown")
  Write (3,*) "Verification results for the TC forecast using dependent datasets"
  Do ifile=1, Nfile
    ! if (ifile<=Nmon) then; Region="VES_N15"; else; Region="VES_S15"; endif
     Region="VES" !"trim(FnameInp(ifile))
     imon=ifile !mod(ifile,Nmon); if (imon==0) imon=Nmon
     Write (3,"(A,A,A,I2.2)") " Input file: ", trim(Region), " Mon: ",imon
     Write (3,"(2(A,F6.1))") " Mean TC-count: ", Mean_Obs(ifile), " Standard Deviation: ", Std_Obs(ifile)
     Write (3,"(A8,20A8)") "ID/Model", (Mdel_nm(imod),imod=2,Nmodels+1)
     Write (3,"(A8,20F8.1)") " ME: ", (ME(ifile,imod), imod=1, Nmodels) 
     Write (3,"(A8,20F8.1)") " MAE: ", (MAE(ifile,imod), imod=1, Nmodels)
     Write (3,"(A8,20F8.1)") " RMSE: ", (RMSE(ifile,imod), imod=1, Nmodels)
     Write (3,"(A8,20F8.1)") " RMAE: ", (RME(ifile,imod), imod=1, Nmodels)
     Write (3,"(A8,20F8.3)") " CORR: ", (Corre(ifile,imod), imod=1, Nmodels)
     Write (3,"(A8,20F8.3)") " Pval: ", (Pvalue(ifile,imod), imod=1, Nmodels)
     Write (3,"(20A8)") ("---------", i=1,6)
  Enddo
End Subroutine Write_Verification_to_TXT
!--------------------------------------------------------------------------------

Subroutine Write_Model_Coeffs_to_NC (Fname)
  Character (Len=*) :: Fname
  Character (Len=255) :: var_ln, t_unit
  Real, Allocatable :: Xlon(:), Xlat(:), Xtime(:), X(:,:,:)
  Integer, Parameter :: NModel=3 ! MLR, LAD, LMV
  Integer :: Mvars

  Mvars=Mbien-1
  Allocate (Xlon(Nfile), Xlat(NModel), Xtime(Mvars+1), X(Nfile,Nmodel,Mvars+1))
  Do ilon=1, Nfile; Xlon(ilon)=ilon; Enddo
  Do ilat=1, Nmodel; Xlat(ilat)=ilat; Enddo
  Do it=1, Mvars+1; Xtime(it)=it-1; Enddo
  Do ifile=1, Nfile
     Do imod=1, Nmodel
        X(ifile,imod,1:Mvars+1) = Beta(ifile,imod,0:Mvars)
     Enddo
  Enddo
  !print*,X
  print*,"Nfile, Nmodel, Mvars=", Nfile, Nmodel, Mvars+1
  t_unit="years since 1961 "
  var_ln="Cac he so cua cac PTHQ b0, b1,... cua cac mo hinh MLR, LAD, LMV"
  Call wr_ncdf_3d (trim(Fname),XLon,XLat,Xtime, Nfile,Nmodel,Mvars+1, X, "Beta", & 
                   trim(t_unit), trim(var_ln),"-","files","Model","heso",Miss)

  Deallocate (Xlon, Xlat, Xtime, X)
End Subroutine Write_Model_Coeffs_to_NC
!--------------------------------------------------------------------------------
Subroutine Write_Verification_Results_to_NC (Fname)
  Character (Len=255) :: var_ln
  Character (Len=*) :: Fname
  Real, Allocatable :: Xlon(:), Xlat(:), X(:,:)
  Integer :: ilon, ilat, it

  Allocate (Xlon(Nfile),Xlat(Nmodels), X(Nfile, Nmodels))
  Do ilon=1, Nfile; Xlon(ilon) = ilon; Enddo      ! So file
  Do ilat=1, Nmodels; Xlat(ilat) = ilat; Enddo    ! 4 Models

  var_ln="Mean Error (ME) forecast of TC-count by different models for Jan to Dec and for North and South VES (Case)"
  X(1:Nfile,1:Nmodels) = ME(1:Nfile,1:Nmodels)
  Call wr_ncdf_2d     (trim(Fname), Xlon, Xlat, Nfile, Nmodels, X, "ME", "Mean Error", "NTC", "Case", "Model", Miss)

  var_ln="Mean Absolute Error (MAE) forecast of TC-count by different models for Jan to Dec and for North and South VES (Case)"
  X(1:Nfile,1:Nmodels) = MAE(1:Nfile,1:Nmodels)
  Call wr_ncdf_2d_apd (trim(Fname), X, "MAE", Nfile, Nmodels, "Mean Absolute Error", "NTC", "Case", "Model", Miss)

  var_ln="Root Mean Square Error (RMSE) forecast of TC-count by different models for Jan to Dec and for North and South VES (Case)"
  X(1:Nfile,1:Nmodels) = RMSE(1:Nfile,1:Nmodels)
  Call wr_ncdf_2d_apd (trim(Fname), X, "RMSE", Nfile, Nmodels, "Root Mean Square Error", "NTC", "Case", "Model", Miss)

  var_ln="Mean Relative Absolute Error (RME) forecast of TC-count by different models for Jan to Dec and for North and South VES (Case)"
  X(1:Nfile,1:Nmodels) = RME(1:Nfile,1:Nmodels)
  Call wr_ncdf_2d_apd (trim(Fname), X, "RME", Nfile, Nmodels, "Mean Relative Absolute Error", "%", "Case", "Model", Miss)

  var_ln="Multi Correlation Coefficients (Corr) of TC-count forecast by different models for Jan to Dec and for North and South VES (Case)"
  X(1:Nfile,1:Nmodels) = Corre(1:Nfile,1:Nmodels)
  Call wr_ncdf_2d_apd (trim(Fname), X, "Corr", Nfile, Nmodels, "Multi Correlation Coefficients", "-", "Case", "Model", Miss)

  var_ln="P-value of Multi Correlation Coefficients (Corr) of TC-count forecast by different models for Jan to Dec and for North and South VES (Case)"
  X(1:Nfile,1:Nmodels) = Pvalue(1:Nfile,1:Nmodels)
  Call wr_ncdf_2d_apd (trim(Fname), X, "Prob", Nfile, Nmodels, "P-value of Multi Correlation Coefficients", "-", "Case", "Model", Miss)
  Deallocate (Xlon, Xlat, X)
End Subroutine Write_Verification_Results_to_NC
!--------------------------------------------------------------------------------
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=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)
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) 
     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
  Real, Allocatable :: X(:) ! X(N*(Mvars+1))
  Real :: S,Q,U,RR,F

  Allocate(X(N*(Mvars+1)))
  X(1:N) = Yorg(1:N)
  Do j=1, Mvars
     i=j*N+1 
     X(i:i+N-1) = Xorg(1:N,j) 
  Enddo
  Call REGRE(X,N,Mvars+1,Beta_tmp, S,Q,U,RR,F)
  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
  Deallocate(X)
End Subroutine Calc_HQTT
!--------------------------------------------------------------------------------
Subroutine Calc_LAD
  Real, Allocatable :: XY(:,:)
  Integer :: i, j, IOPT

  Allocate(XY(1:Nmax,1:Mmax+5))
  Do j=1, Mvars
     XY(1:N,j)=Xorg(1:N,j)
  Enddo
  XY(1:N,Mvars+1)=Yorg(1:N)
  IOPT=0; Call RLLAV (XY,N,Mvars,IOPT, Beta_tmp)
  Do i=1, N
     Yhq(i) = 0
     Do j=1, Mvars
        Yhq(i) = Yhq(i)+Beta_tmp(j)*Xorg(i,j)  ! Beta_tmp(1..Mvars)=A1,...Am
     Enddo
     Yhq(i) = Yhq(i)+Beta_tmp(Mvars+1)         ! Beta_tmp(Mvars+1)=A0
  Enddo
  Deallocate(XY)
End Subroutine Calc_LAD
!--------------------------------------------------------------------------------
Subroutine Calc_LMV
  Real, Allocatable :: XY(:,:)
  Integer :: i, j, IOPT

  Allocate(XY(1:Nmax,1:Mmax+5))
  Do j=1, Mvars
     XY(1:N,j)=Xorg(1:N,j)
  Enddo
  XY(1:N,Mvars+1)=Yorg(1:N)
  IOPT=0; Call RLLMV (XY,N,Mvars,IOPT, Beta_tmp)
  Do i=1, N
     Yhq(i) = 0
     Do j=1, Mvars
        Yhq(i) = Yhq(i)+Beta_tmp(j)*Xorg(i,j)  ! Beta_tmp(1..Mvars)=A1,...Am
     Enddo
     Yhq(i) = Yhq(i)+Beta_tmp(Mvars+1)         ! Beta_tmp(Mvars+1)=A0
  Enddo
  Deallocate(XY)
End Subroutine Calc_LMV
!--------------------------------------------------------------------------------
Subroutine Calc_Verification (Yorg,Yhq,N, ME, MAE, RMSE, RME, Corr, Pvalue)
! Input: Yorg(1:N) is an array stored the origin values of Y
!        Xorg(1:N,1:M) is an array stored the origin values of X
!        Beta(0:M) is an array stored the regression coeficients B0, B1,..., BM
! Output: ME is Mean Error
!         MAE is Mean Absolute Error
!         RMSE is Root Mean Square Error
!         RME is Relative Mean Error
!         Corr, Pvalue are Correlation Coefficient and P-value (t-test)
Integer :: N, i, j, k, Nc, Nc0
Real    :: Yorg(N), Yhq(N), Error
Real    :: ME, MAE, RMSE, RME, Corr, Pvalue

  ME=0; MAE=0; RMSE=0; RME=0; Nc=0; Nc0=0 
  Do i=1, N
     if (Yhq(i)/=Miss.AND.Yorg(i)/=Miss) then
        Nc=Nc+1
        Error= Yhq(i)-Yorg(i)
        ME   = ME   + Error
        MAE  = MAE  + ABS(Error)
        RMSE = RMSE + Error*Error
        if (Yorg(i)/=0) then; Nc0=Nc0+1; RME  = RME  + ABS(Error/Yorg(i)); endif
     endif
  Enddo
  ME = ME/Real(Nc); MAE = MAE/Real(Nc); RME = RME/Real(Nc0)*100.; RMSE = Sqrt(RMSE/Real(Nc))
  Call Correlation (Yorg, Yhq, N, Corr, Pvalue, Miss)
End Subroutine Calc_Verification
!--------------------------------------------------------------------------------
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_neural.f90"
Include "SUB/Stdev.f90"
Include "SUB/Ave.f90"
!--------------------------------------------------------------------------------
END
Include "SUB/Student_Prob.f90"
!-------------------------------------------------------------------------------------------