Use NetCDF
   Integer, Parameter :: Yr1=1982, Yr2=2011, Nmon=12 
   Integer, Parameter :: NDr_lev=4, Nlen=3
   Character (Len=3),Parameter :: Dr_lev(NDr_lev)=(/"mil","mod","sev","ext"/)
   Real, parameter    :: Miss=-99.0
   Character (Len=255) :: Fname, Inp_Dir, Mask_File, Out_Dir

   Real, Allocatable, Dimension (:,:,:)   :: PC, POD, TS, HSS ! (Nlon,Nlat,Nlead)
   Real, Allocatable, Dimension (:,:,:,:,:) :: beta             ! (Nlon,Nlat,Nlead,Ngroup,NVAR+1)
   INTEGER, ALLOCATABLE :: NG(:)                              ! Dung luong mau thuc te
   REAL, ALLOCATABLE :: X(:,:), X1(:,:), MEAN(:,:), C(:,:), C1(:,:), COEF(:,:)
   REAL, ALLOCATABLE :: F(:,:)
   INTEGER, ALLOCATABLE :: LG(:),LFSC(:)

   Real, Allocatable :: Regs(:,:), Lon(:), Lat(:)             ! Lon, lat and regions

   Integer :: Nlon,Nlat,Nlead, ireg,ilon,ilat,ist, Nt,it, icnt
   Integer :: Nmods,Ncase,ilen,idr
   Integer :: Ntime,M, i,j,k, ilead, ndf, ier, IO
   Real    :: aa,bb,cc,dd

   INTEGER :: NGROUP, NVAR, NMAX, N
   Integer :: icount(6), ic1(6)
   Logical :: succes

   !------------------------------------------------------------------------------
   Ntime=(Yr2-Yr1+1)*Nmon   ! Maximum timeseries
   Nlead=5                  ! All leadtimes for baseline period
   Mask_File    = "/home/tanpv/Data/AWC/Mask_VN_0.25_Corrected.nc"
   Call rd_ncdf_2d (Mask_File,Lon,Lat,Nlon,Nlat,Regs,"reg")

   print*,Nlon,Nlat,Nlead

   Allocate (PC(Nlon,Nlat,Nlead),POD(Nlon,Nlat,Nlead),TS(Nlon,Nlat,Nlead),HSS(Nlon,Nlat,Nlead))
   PC=Miss; POD=Miss; TS=Miss; HSS=Miss;
   ilen=1            ! Drought timescale: 1: 1-month, 3: 3-month, 6: 6-month
   idr =2            ! Drought level:     1: Moderate, 2: Severe, 3:Extreme
   Out_Dir = "/home/tanpv/Data/562_VNU/0.25_VN/4Discr/"

   write (Fname,"(A,I2.2,A)") trim(Out_Dir)//"Dr_evt_and_model_out_4_discr_logis_Dr_"//trim(Dr_lev(idr))//"_len_",ilen,".dat"
   print*,trim(Fname)
   Open (1,file=Fname,status="Old",form='unformatted')

   IO=0; icnt=0; icount(:)=0; ic1(:)=0
      Do while (io==0)
         !if (IO/=0) exit
         Read (1,IOSTAT=io) ilead, ireg !ilon, ilat
print*,ilead, ireg !ilon, ilat
         if (IO/=0) exit
         Read (1) Nt, Nvar, Ngroup  ! Sample size, Num of var, Num of group 
         M=Nvar                     ! 
         icnt=icnt+1
         ii=mod(ilead,Nlead); if (ii==0) ii=Nlead
         icount(ii)=icount(ii)+1
!print*,icount
         if (allocated (LG)) Deallocate (LG); if (allocated (X))  Deallocate (X)
         if (allocated (NG)) Deallocate (NG); if (allocated (MEAN)) Deallocate (MEAN)
         if (allocated (COEF)) Deallocate (COEF); if (allocated (C)) Deallocate (C)
         if (allocated (F)) Deallocate (F); if (allocated (LFSC)) Deallocate (LFSC)

         Allocate (NG(Ngroup),LG(Nt), X(Nt,M))
         Allocate (COEF(Ngroup,M+1))
         Allocate (F(Nt,NGROUP), LFSC(Nt))
                  
         if (icnt==1) then; Allocate (beta(Nlon,Nlat,Nlead,Ngroup,M+1)); beta=Miss; endif

         Read (1) NG
         Read (1) LG
         Read (1) X

         Call Discr_2023(NGROUP,Nt,NVAR,NG,X,COEF,LFSC)
         !Do k=1,NGROUP
         !   print"(I3,10F8.3)",k, (COEF(k,j),j=1,Nvar+1)
         !Enddo
         ! Handle errors:
         succes=.true.
         Do k=1,NGROUP
            Do j=1,Nvar+1
               if (IsNaN(COEF(k,j))) then
                  succes=.False.; exit
               endif
            Enddo
         Enddo
         if (succes) then
            ic1(ii)=ic1(ii)+1
            ! Saving COEF to write out:
            Do k=1,Ngroup
               beta(ilon,ilat,ilead,k,1:M+1)=COEF(k,1:M+1)
            Enddo
            ! Tinh ham phan biet cho tap mau phu thuoc
            Do I=1,Nt ! 1
               Do K=1,NGROUP
                  F(I,K) = COEF(K,1)
                  Do J=1,NVAR
                     F(I,K) = F(I,K) + COEF(K,J+1)*X(I,J)
                  Enddo
               Enddo
               LFSC(I) = 1
               Temp = F(I,1)
               Do K=2,NGROUP
                  if (F(I,K)>Temp) THEN
                     Temp = F(I,K)
                     LFSC(I) = K
                  endif
               Enddo
            Enddo
            !Do I=1,Nt ! 1
            !   print '(2F20.4,2I4)', (F(I,K),K=1,NGROUP), LFSC(I)
            !Enddo

            aa=0; bb=0;cc=0;dd=0
            Do i=1,Nt
               if (LFSC(i)==1.AND.LG(i)==1) aa=aa+1
               if (LFSC(i)==1.AND.LG(i)==0) bb=bb+1
               if (LFSC(i)==0.AND.LG(i)==1) cc=cc+1
               if (LFSC(i)==0.AND.LG(i)==0) dd=dd+1
            Enddo
            print*,"aa,bb,cc,dd: ",aa,bb,cc,dd
            print*,(aa+dd)/(aa+bb+cc+dd),aa/(aa+cc); 
            pc=aa+bb+cc+dd; pc=(aa+dd)/pc
            PC(ilon,ilat,ilead) =(aa+dd)/(aa+bb+cc+dd)
            POD(ilon,ilat,ilead)=aa/(aa+cc);
            TS(ilon,ilat,ilead) =aa/(aa+bb+cc);
            HSS(ilon,ilat,ilead)=2.0*(aa*dd-bb*cc)/((aa+cc)*(cc+dd)+(aa+bb)*(bb+dd))
         else
            !print*, "Not success for ",ilead,ilon,ilat
         endif
      Enddo
      print*,icount(:)
      print*,ic1(:)
   CONTAINS
    !------------------------------------------------------------------------------------

   Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
   Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
   Include "/home/tanpv/Mylibs/Sub_My_Discr_2023.f90"

   END