Use NetCDF

      Real, Parameter :: Miss = -99.0
      Real, allocatable :: xlon(:), xlat(:), xtime(:)
      Real, Allocatable :: Xmean(:,:), SS(:,:)
      Integer :: Nlon, Nlat, Ntime
      Character (Len=100) :: Fname, Fname1
      Logical :: Status
!------------------------------------------------------------------------------------
      Real, Parameter :: Yr1=1976, Yr2=2005, Nyear=Yr2-Yr1+1, Nmon=12
      Real, allocatable :: PDSI(:,:,:), SPI12(:,:,:), SPEI12(:,:,:), Mask(:,:)
      Real, allocatable :: All_IDX(:,:,:,:)
      Real, allocatable :: Dr_IDX(:,:,:)  ! (Nst,Ntime,Ndr_IDX)
      Integer :: Nst, iidx
      Integer, Parameter :: NDr_idx=3 ! PDSI, SPEI12, SPI12
      Character (Len=5), Parameter :: IDX_nm(NDr_idx)= &
                                   (/"PDSI","SPI12","SPE12"/)
      Real, Parameter :: Thre(NDr_idx) = (/-3., -1.5, -1.5/)
      Real, Allocatable :: Ave_DrMon(:,:,:), Ave_DrSev(:,:,:)  ! So thang han va Severity (Nlon, Nlat, NDr_idx)
      Real, Allocatable :: Geo(:,:)  ! Geo Ext (Ntime, NDr_idx)
!------------------------------------------------------------------------------------
      Fname="/home/tanpv/Progs/Dr_SEA/Results/Dr_PDSI_SPEI_SPI.nc"
      Call Read_Drough_IDX (Fname)
      Allocate (Ave_DrMon(Nlon, Nlat, NDr_idx), Ave_DrSev(Nlon, Nlat, NDr_idx), Geo(Ntime, NDr_idx))

      Call Calc_Characteristics
   
      Fname="Dr_Chars.nc"; Fname1="Dr_Geo.nc"
      Call Write_to_NC (Fname,Fname1)

      CONTAINS
!------------------------------------------------------------------------------------
      Subroutine Write_to_NC (Fname1,Fname2)
         Character (Len=*) :: Fname1,Fname2
         Character (Len=255) :: var_ln, t_unit
         Real, Allocatable :: X(:,:), Y(:)
         Integer :: iidx

         Allocate (X(Nlon,Nlat))
         Do iidx=1, NDr_idx
            X(:,:) = Ave_DrMon(:,:,iidx)
            var_ln="Mean of Number of Drought Months by "//trim(IDX_nm(iidx))
            if (iidx==1) then
               Call wr_ncdf_2d (Fname1,Xlon,xlat, Nlon, Nlat, X, "DMo_"//trim(IDX_nm(iidx)), &
                     trim(var_ln), "-", "lon", "lat", Miss)
            else
               Call wr_ncdf_2d_apd (Fname1, X, "DMo_"//trim(IDX_nm(iidx)), Nlon, Nlat, &
                     trim(var_ln), "-", "lon", "lat", Miss)
            endif
         Enddo
         Do iidx=1, NDr_idx
            X(:,:) = Ave_DrSev(:,:,iidx)
            var_ln="Mean of Drought Severity by "//trim(IDX_nm(iidx))
            Call wr_ncdf_2d_apd (Fname1, X, "DSev_"//trim(IDX_nm(iidx)), Nlon, Nlat, &
                     trim(var_ln), "-", "lon", "lat", Miss)
         Enddo
         Deallocate (X);  Allocate (Y(NDr_idx))
         !Do iidx=1, NDr_idx
         !   Y(:) = Geo(:,iidx)
            var_ln="Drought Geographic Extent by different drought index" !//trim(IDX_nm(iidx))
         !   Call wr_ncdf_1d_apd (Fname1, X, "Geo_"//trim(IDX_nm(iidx)), Ntime, &
         !          trim(var_ln), "%", "mon", Miss )
         !Enddo
         Y(1)=1; Y(2)=2; Y(3)=3
         Call wr_ncdf_2d (Fname2, xtime, Y, Ntime, NDr_idx, Geo, "Geo", trim(var_ln), "%", "mon", "IDX", Miss)
      End Subroutine Write_to_NC
!------------------------------------------------------------------------------------
      Subroutine Calc_Characteristics
         Integer :: ilon,ilat, ist, it, iidx
         Real :: idro, isev, igeo, Nlandpoint
         Ave_DrMon=Miss; Ave_DrSev=Miss; Geo=Miss
         Do iidx=1, NDr_idx
            Nlandpoint=0
            Do ilon=1, Nlon
               Do ilat=1, Nlat
                  if (Mask(ilon,ilat)>0) then
                     idro=0; isev=0; Nlandpoint=Nlandpoint+1
                     Do it=1, Ntime
                        if (All_IDX(ilon,ilat,it,iidx)<=Thre(iidx)) then
                           idro = idro + 1  ! So thang han
                           isev = isev + All_IDX(ilon,ilat,it,iidx) ! Severity 
                        endif
                     Enddo
                     Ave_DrMon(ilon,ilat,iidx) = idro/Real(Nyear)   ! Trung binh so thang han
                     Ave_DrSev(ilon,ilat,iidx) = isev/Real(Nyear)   ! Trung binh severity
                  endif
               Enddo
            Enddo
            ! Geo:
            Do it=1, Ntime
               igeo=0
               Do ilon=1, Nlon
                  Do ilat=1, Nlat
                     if (Mask(ilon,ilat)>0.AND.All_IDX(ilon,ilat,it,iidx)<=Thre(iidx).AND.All_IDX(ilon,ilat,it,iidx)/=Miss) igeo=igeo+1
                  Enddo
               Enddo
               Geo(it,iidx) = igeo/Nlandpoint*100
            Enddo
         Enddo
      End Subroutine Calc_Characteristics
!------------------------------------------------------------------------------------
      Subroutine Read_Drough_IDX (Fname)
         Integer :: ist, ilon, ilat, it, Nt, imon, imon1, iyr
         !Integer :: Nyear, Nmon
         Character (Len=*) :: Fname
         Real :: tmp, xmean(12), ic(12)
         !PDSI(:,:,:), SPI12(:,:,:), SPEI12(:,:,:), Mask(:,:)
         !xlon(:), xlat(:), xtime(:), Nlon, Nlat, Ntime

         Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, PDSI, "ori_pdsi")
         Deallocate (xlon, xlat, xtime)
         Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, SPEI12, "SPEI_12")
         Deallocate (xlon, xlat, xtime)
         Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, SPI12, "SPI_12")
         Deallocate (xlon, xlat)
         Call rd_ncdf_2d (trim(Fname), xlon,xlat, nlon,nlat,Mask,"land")
         Ntime = Nt-Nmon
         Allocate (All_IDX(Nlon,Nlat,Ntime,NDr_idx))
         All_IDX(:,:,1:Ntime,1) = PDSI(:,:,Nmon+1:Nt) 
         All_IDX(:,:,1:Ntime,2) = SPI12(:,:,Nmon+1:Nt) 
         All_IDX(:,:,1:Ntime,3) = SPEI12(:,:,Nmon+1:Nt) 
  
         ist=0
         Do ilat=1, Nlat
            Do ilon=1, Nlon
               if (Mask(ilon,ilat)>0) ist=ist+1
            Enddo
         Enddo
         Nst=ist
         Allocate (Dr_IDX(Nst,Ntime,NDr_idx))
         ist=0
         Do ilat=1, Nlat
            Do ilon=1, Nlon
               if (Mask(ilon,ilat)>0) then
                  ist=ist+1
                  Dr_IDX(ist,1:Ntime,1)=All_IDX(ilon,ilat,1:Ntime,1)
                  Dr_IDX(ist,1:Ntime,2)=All_IDX(ilon,ilat,1:Ntime,2)
                  Dr_IDX(ist,1:Ntime,3)=All_IDX(ilon,ilat,1:Ntime,3)
               endif
            Enddo
         Enddo
      End Subroutine Read_Drough_IDX
!------------------------------------------------------------------------------------
      Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
      Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
      Include "/home/tanpv/Mylibs/Ave.f90"
!------------------------------------------------------------------------------------
      END