Use NetCDF

   Real, Parameter :: Miss=-99.0
   Integer, Parameter :: IMiss=-99, Nvar=5, N365=365 
   Integer, Parameter :: Mwd=15           ! Running mean window in days
   Integer, Parameter :: Swd=7            ! Change window in day
   
   Character (Len=5), Parameter :: vname(Nvar)=(/"pev","e","swvl1","swvl2","swvl3"/)
   Character (Len=25), Parameter :: var_ln(Nvar)=(/"Potential evaporation (m)","lwe_thickness_of_water_evaporation_amount (m)", &
                                    "Volumetric soil water layer 1 (m**3 m**-3)","Volumetric soil water layer 2 (m**3 m**-3)", &
                                    "Volumetric soil water layer 3 (m**3 m**-3)"/)
   Real, Allocatable, Dimension (:)     :: Lon, Lat, time
   Integer                              :: Nlon, Nlat, Ntime
   Real, Allocatable, Dimension (:,:,:) :: PET, ET, sw1, sw2, sw3 
   Real, Allocatable, Dimension (:,:,:) :: Pet_wd, Et_wd ! , Pet_31, Et_31
   Real, Allocatable, Dimension (:,:,:) :: ESI_wd        ! , ESI_31
   Real, Allocatable, Dimension (:,:,:) :: dESI_wd       ! Change in ESI over Swd window steps
   Real, Allocatable, Dimension (:,:,:) :: zESI          ! standardized ESI change anomalies
   Real, Allocatable, Dimension (:,:,:) :: fd_evt        ! FD event (0/1)

   Character (Len=255) :: OutDir, Fname, t_unit, DirIn, Str, time_unit, Fname1, Fname2
   Character (Len=1) :: Ch
   Real*8 :: scale_factor, add_offset, missing_value, fill_value 
   Integer :: No, iday, ievt
   !----------------------------------------------------------------------------------------------

   Fname="flash_drough_26-04.csv"
   Open (1,file=Fname,status="old")
   Nlon=26; Nlat=31; Ntime=7300
   Allocate (Lon(Nlon), Lat(Nlat), time(Ntime))
   Allocate (ESI_wd(Nlon,Nlat,Ntime),zESI(Nlon,Nlat,Ntime))
   Allocate (fd_evt(Nlon,Nlat,Ntime))
   Read (1,*)
   Do ilat=1,Nlat !ilon=1,Nlon
      Do ilon=1,Nlon !ilat=1,Nlat
         Do it=1,Ntime
            Read (1,*) No, iday, Lon(ilon),Lat(ilat),ESI_wd(ilon,ilat,it),zESI(ilon,ilat,it),fd_evt(ilon,ilat,it)
         Enddo
         !print*,Lon(ilon),Lat(ilat)
      Enddo
!pause
   Enddo
   print*,"Reading ... Done"
   Fname="esi_zesi_fd_soctrang_2001-2020.nc"
   Call Write_data_to_file (Fname)
   CONTAINS

   !----------------------------------------------------------------------------------------------
   Subroutine Write_data_to_file (Fname) 
      Character (Len=*) :: Fname
      Integer :: it,ifile,ireg,Ndays_60y
      Real, Allocatable :: X3d(:,:,:),xday(:),X2d(:,:),xreg(:)
      Character (Len=150) :: var_ln, t_unit, lon_nm, lat_nm, z_nm, t_nm
      Character (Len=10) :: var_unit

      Do it=1,Ntime; time(it)=it-1; Enddo

      var_unit="-"; t_unit="days since 2001-01-01 ";  lon_nm="lon"; lat_nm="lat"; t_nm="time"
      write(var_ln,"(A,I2.2,A)") "ESI with ",Mwd,"-day window running mean of PET and ET"
      Call wr_ncdf_3d (Fname,Lon,Lat,time,Nlon,Nlat,Ntime,ESI_wd,"esi",t_unit,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
      write(var_ln,"(A,I2.2,A)") "Standardized ESI change anomalies (zESI) over ",Swd,"-day window intervals"
      Call wr_ncdf_3d_apd (Fname,zESI,"zesi",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
      write(var_ln,"(A,I2.2,A)") "FD event: 0/1"
      Call wr_ncdf_3d_apd (Fname,fd_evt,"fd",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
   End Subroutine Write_data_to_file
   !----------------------------------------------------------------------------------------------

   Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
   Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
   Include "/home/tanpv/Mylibs/Ave.f90"
   Include "/home/tanpv/Mylibs/Stdev.f90"
   Include "/home/tanpv/Mylibs/Xmax.f90"
   Include "/home/tanpv/Mylibs/Xmin.f90"
   Include "/home/tanpv/Mylibs/Xsum.f90"
   Include "/home/tanpv/Mylibs/Sub_Percentile.f90"
   Include "/home/tanpv/Mylibs/Percentile_rank.f90"
   END