Use NetCDF
      Real, Parameter :: Miss=-99.0
      Integer, Parameter :: Nidx=2, NEnso=4, Yr1=1981, Yr2=2020, 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"
      Fname="/home/tanpv/Data/SEACLID/ClimIDX_2022/Clim_IDX_2022.nc"
      Call Read_ClimIDX (Fname)
      Allocate (ENSO(Ntime,NEnso),Xlon(Ntime),Xlat(NEnso))
      Fname="/home/tanpv/Progs/Dr_4Books/Res/ONI_ENSO_1981-2020.nc"
      !Call Write_ENSO_timeseries (Fname) 
      Call Write_ENSO (Fname)

      CONTAINS
      !--------------------------------------------------------------------------------------
      Subroutine Write_ENSO_timeseries (Fname)
         Integer :: ilon, ilat, it, ist
         Character (Len=*) :: Fname
         Character (Len=255) :: t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm
         Real, Allocatable :: ENSO_out(:,:,:) 
         Real :: Dummy(1)
         !----------------------------------------------------------------------
         ! write ONI as Enso definition: ONI<=0.5: La Nina; >=0.5: El Nino; in between this interval is Neutral
         Allocate (ENSO_out(1,NEnso,Ntime))
         ENSO_out(:,:,:)=0.0
         Do ilon=1, Ntime; Xlon(ilon)=ilon; Enddo        ! Mang time
         Xlat(1)=-1; Xlat(2)=0; Xlat(3)=1; Xlat(4)=3     ! All=3; El Nino=1; La Nina=-1; Neutral=0
         Dummy(1)=1.
         lon_nm="none"; lat_nm="ENSO_phase"; t_nm="month"
         t_unit="months since 1981-01-01"; var_unit="C"
         var_ln="ONI used as ENSO classification"
         Do it=1, Ntime
            ENSO_out(1,4,it)=ClimIDX(1,it)           ! All
            if (ClimIDX(1,it)<=-0.5) then      ! La Nina
               ENSO_out(1,1,it)=ClimIDX(1,it)
            else if (ClimIDX(1,it)>=0.5) then  ! El Nino
               ENSO_out(1,3,it)=ClimIDX(1,it)
            else
               ENSO_out(1,2,it)=ClimIDX(1,it)        ! Neutral
            endif
         Enddo
         Call wr_ncdf_3d (trim(Fname),Dummy,Xlat,Xlon,1,4,Ntime,ENSO_out,"ONI", &
                    t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      End Subroutine Write_ENSO_timeseries
      !--------------------------------------------------------------------------------------
      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