Use NetCDF
    Use neural

    type(network) :: testNet

    Character (Len=255) :: Fname

    Integer :: ireg, ilead, Nreg, Nlead
    Integer :: Nlayer                             ! Num of layers including Input, output and hidden laysers
    Integer :: Nepoch                             ! Num of epochs
    Integer :: Nnode_max                          ! Maximum number of hidden nodes
    double precision :: Learning_Rate
    Integer, Allocatable :: Layer_size(:)                      ! (Nlayer)

    double precision, Allocatable :: tmp(:)
    Real :: a,b,c,d,PC,POD, Thre, PC0,POD0

    Integer :: ilr, inode, itest, Ntest

    !-----------------------------------------------------------------------------
    Fname="/home/tanpv/Data/562_VNU/0.25_VN/Dr_IDXs_Ope/Reg_ens_Dr_IDX_1983-2020_spei_01.nc"
    Call Read_data

    Nepoch        = 10000
    Learning_Rate = 0.01d00
    Nlayer        = 5
    Nnode_max     = 15
    Allocate (Layer_size(Nlayer))
    Layer_size(1) = Mxvar; Layer_size(Nlayer) = Myvar
 
    Allocate (X_inp(Nspl,Mxvar),Y_out(Nspl,Myvar),Yest(Nspl,Myvar))
    ilead = 1
    Do ireg=1,Nreg
       Y_out(:,1) = Obs (ireg,ilead,:)
       X_inp(:,1) = Ens (ireg,ilead,:)
       X_inp(:,2) = Ens1(ireg,ilead,:)
       PC0=0.0; POD0=0.0
       Do inode=Nnode_max,1,-1
          Do ilr=2,Nlayer-1
             Layer_size(ilr) = inode
          Enddo
          call neural_init_network(testNet, Layer_size) ! (/ 2, 15, 15, 1 /))

          call neural_train_network_file_TP(testNet, Fname, Nepoch, Learning_Rate)
          !call neural_save_network(testNet, "network.dat")
          !call neural_load_network(testNet, "network.dat")

          if (Allocated (X)) Deallocate (X); if (Allocated (tmp)) Deallocate (tmp)
          Allocate (X(Mxvar),tmp(Myvar))
          a=0.;b=0.;c=0.;d=0.; Thre=-0.5
          Do i=1,Nspl
             X(1:Mxvar)=X_inp(i,1:Mxvar)
             tmp(:) = neural_use_network(testNet,X)
             Do j=1,Myvar
                tmp(j) = tmp(j)*(Ymax(j)-Ymin(j))+Ymin(j)
                Y_out(i,j) = Y_out(i,j)*(Ymax(j)-Ymin(j))+Ymin(j)
             Enddo
             print "(10F7.3)",tmp(:),Y_out(i,:) !,Obs (1,1,i)
             Where (tmp<=Thre)
                tmp = 1.0
             Elsewhere
                tmp = 0.0
             Endwhere
             Where (Y_out(i,:)<=Thre)
                Y_out(i,:) = 1.0
             Elsewhere
                Y_out(i,:) = 0.0
             Endwhere
             if (tmp(1) == 1. .AND. Y_out(i,1)==1.) a=a+1.   ! Obs=1; Fcs=1
             if (tmp(1) == 0. .AND. Y_out(i,1)==0.) d=d+1.   ! Obs=0; Fcs=0
             if (tmp(1) == 0. .AND. Y_out(i,1)==1.) c=c+1.   ! Obs=1; Fcs=0
             if (tmp(1) == 1. .AND. Y_out(i,1)==0.) b=b+1.   ! Obs=0; Fcs=1
          Enddo
          PC  = (a+d)/(a+b+c+d)
          POD = (a+c)/(a+b+c+d)
          print*,a,b,c,d
          print*,PC,POD
          if (POD>POD0) then
             POD0=POD; PC0=PC; itest=inode
          endif
          Y_out(:,1) = Obs (ireg,ilead,:)
          X_inp(:,1) = Ens (ireg,ilead,:)
          X_inp(:,2) = Ens1(ireg,ilead,:)

       Enddo
       print*,"The best: ",itest,PC0,POD0; pause
    Enddo 
CONTAINS

Subroutine Read_data
!    integer :: ierr, input_size, output_size, num_lines, i, j, num_data
!    Integer :: ireg, ilead, Nreg, Nlead, Ntime
    Real, allocatable :: xreg(:), xlead(:), xtime(:) !, Obs(:,:,:), Ens(:,:,:), Ens1(:,:,:)

    Call rd_ncdf_3d (Fname,xreg,xlead,xtime,Nreg, Nlead, Ntime,Obs,"obs")
    Deallocate (xreg,xlead,xtime)
    Call rd_ncdf_3d (Fname,xreg,xlead,xtime,Nreg, Nlead, Ntime,Ens,"ens")
    Deallocate (xreg,xlead,xtime)
    Call rd_ncdf_3d (Fname,xreg,xlead,xtime,Nreg, Nlead, Ntime,Ens1,"ens1")
    Deallocate (xreg,xlead,xtime)

    Mxvar = 2; 
    Myvar = 1; 
    Nspl  = Ntime
End Subroutine Read_Data
Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
End