Use NetCDF
      Real, Parameter :: Miss=-99.0
      Integer, Parameter :: Nidx=2, NEnso=4, Yr1=1960, Yr2=2019, Nmon=12
      Character (Len=255) :: Fname
      Character (Len=10) :: Clim_nm(Nidx)=(/"oni","soi"/)
      Real, Allocatable :: ClimIDX(:,:), ENSO(:,:), Xlon(:), Xlat(:)
      Integer :: Ntime 
      !--------------------------------------------------------------------------------------
      Fname="/work/projects/G09/SEACLID/ClimIDX/Clim_IDX_2020.nc"
      Call Read_ClimIDX (Fname)
      Allocate (ENSO(Ntime,NEnso),Xlon(Ntime),Xlat(NEnso))
      Fname="ONI_ENSO_1960-2019.nc"
      Call Write_ENSO (Fname)

      CONTAINS
      !--------------------------------------------------------------------------------------
      Subroutine Write_ENSO (Fname)
         Integer :: ilon, ilat, it, ist
         !Integer, Parameter :: NEnso=4
         !Real, Dimension (N_ClimIDX, Ntime) :: ClimIDX

         !Real :: ENSO(Ntime,NEnso), Xlon(Ntime), Xlat(NEnso)
         Character (Len=*) :: Fname
         !----------------------------------------------------------------------
         ! write ONI as Enso definition: ONI<=0.5: La Nina; >=0.5: El Nino; in between this interval is Neutral
         ENSO(:,:)=0.0
         Do ilon=1, Ntime; Xlon(ilon)=ilon; Enddo        ! Mang Lon
         Xlat(1)=-1; Xlat(2)=0; Xlat(3)=1; Xlat(4)=3     ! All=3; El Nino=1; La Nina=-1; Neutral=0
         Do it=1, Ntime
            ENSO(it,4)=ClimIDX(1,it)           ! All
            if (ClimIDX(1,it)<=-0.5) then      ! La Nina
               ENSO(it,1)=ClimIDX(1,it)
            else if (ClimIDX(1,it)>=0.5) then  ! El Nino
               ENSO(it,3)=ClimIDX(1,it)
            else
               ENSO(it,2)=ClimIDX(1,it)        ! Neutral
            endif
         Enddo

         Call wr_ncdf_2d (trim(Fname), Xlon, Xlat, Ntime, NEnso, ENSO, &
             "ONI", "ONI used as ENSO classification", "C", "Mon","Enso",Miss)

      End Subroutine Write_ENSO
      !--------------------------------------------------------------------------------------
      Subroutine Read_ClimIDX (Fname)
       Character (Len=*) :: Fname
       Real, Allocatable :: Year(:),Month(:), tmp(:,:)
       Integer :: Nyr, Nmn, iidx, iyr, i,j, k

       Ntime=(Yr2-Yr1+1)*Nmon
       Allocate (ClimIDX(Nidx,Ntime)) ! ,Nidx))

       !print*,Nidx,Ntime,Yr1,Yr11,Yr2; pause
       Do iidx=1,Nidx
          Call rd_ncdf_2d (Fname,Year,Month,Nyr,Nmn,tmp, trim(Clim_nm(iidx)))
          j=0
          Do iyr=Yr1,Yr2
             Do i=1,Nyr
                if (Year(i)==iyr) then
                   Do k=1,Nmon
                      j=j+1; ClimIDX(iidx,j)=tmp(i,k)
                   Enddo
                endif
             Enddo
          Enddo
          Deallocate (Year,Month,tmp)
          !print "(12F7.3)",Y(iidx,:)
          !print*,Clim_nm(iidx);pause
       Enddo
      End Subroutine Read_ClimIDX
      !--------------------------------------------------------------------
      Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
      Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"

      END