Use NetCDF

   Real, Parameter :: Miss=-99.0
   Integer, Parameter :: IMiss=-99, Nfiles=9, Nreg=7, Pentad=5, N365=365, Nwd=5  ! Nwd=Window data
   Integer, Allocatable :: time(:)

   Integer :: i, LenSt, sys
   Character (Len=255) :: OutDir, Fname, t_unit, DirIn, Str, time_unit, Fname1, Fname2
   Character (Len=255), Allocatable :: Files(:) 
   Character (Len=12) :: date0, timeres, idate0, idate1
   Character (Len=1) :: Ch
   Integer :: yy1,mm1,dd1,hh1,yy2,mm2,dd2,hh2
   Real*8 :: factor, offset, Missing, fill_value
   Real, Allocatable :: Prec(:,:,:), lon(:), lat(:), lev(:), xtime(:), Regs(:,:)
   Real, Allocatable, Dimension (:,:,:) :: sw1,sw2,sw3,sw01,sw02,sw03   ! (Nlon,Nlat,Ntime) 
   Real, Allocatable, Dimension (:,:)   :: sw1_Reg,sw2_Reg,sw3_Reg,sw01_Reg,sw02_Reg,sw03_Reg ! (Nreg,Ntime)
   
   Real, Allocatable, Dimension (:,:,:) :: psw1,psw2,psw3,psw01,psw02,psw03   ! (Nlon,Nlat,Ntime) 
   Real, Allocatable, Dimension (:,:)   :: psw1_Reg,psw2_Reg,psw3_Reg,psw01_Reg,psw02_Reg,psw03_Reg ! (Nreg,Ntime)
   !----------------------------------------------------------------------------------------------
   ! Precipitation (percentile), T2m, Tx, Tm, Evaporation:
   Real, Allocatable, Dimension (:,:,:) :: TP, TPp, T2, Tx, Tm, Eva      ! (Nlon,Nlat,Ntime)
   Real, Allocatable, Dimension (:,:)   :: TP_Reg, TPp_Reg, T2_Reg, Tx_Reg, Tm_Reg, Eva_Reg   ! (Nreg,Ntime)
   !----------------------------------------------------------------------------------------------
   Real, Allocatable, Dimension (:,:,:,:) :: X4d
   Integer, Allocatable, Dimension (:) :: date_m1   ! (Nfcs)
   Integer :: Nlon, Nlat, Nlev, Ntime, Ndays1, Ndays, Npentad, Nlon0, Nlat0    ! Nfcs = So lan du bao, date_fcs = yyyymmdd cua ngay DB

   Integer :: imem, ifcs, iy, im, id

   im=0; imem=365; dd1=21900; k=0; ifcs=0
   do i=1,imem
      k=0
      Do it=i,dd1,imem
         do j=it-2,it+2
            if (j>0.AND.j<dd1) then
               k=k+1
    !           print*,"i,k,it,j=",i,k,it,j
            endif
         enddo
      Enddo
      !pause
   enddo
   !print*,"ifcs=",ifcs
!pause

   Call Read_Regs ("/home/tanpv/Data/ERA5/Mask_VN_0.25.nc") !VN.grid_0.25.nc")

   Call Read_sw_data
   Call Calculating_SM_Percentile
   OutDir="/home/tanpv/Progs/Flash_Dr/NC/"
   Call Write_data_to_file (trim(OutDir)//"percentile_sw_1961-2020_daily.nc",trim(OutDir)//"percentile_sw_1961-2020_daily_Reg.nc")
   CONTAINS
   !----------------------------------------------------------------------------------------------
   Subroutine Calculating_SM_Percentile
      Integer :: ilon, ilat, it, ireg
      Real :: X(Ntime), Y(Ntime)

      Allocate (psw1(Nlon,Nlat,Ntime),psw2(Nlon,Nlat,Ntime),psw3(Nlon,Nlat,Ntime))
      Allocate (psw01(Nlon,Nlat,Ntime),psw02(Nlon,Nlat,Ntime),psw03(Nlon,Nlat,Ntime))
      Allocate (psw1_Reg(Nreg,Ntime),psw2_Reg(Nreg,Ntime),psw3_Reg(Nreg,Ntime))
      Allocate (psw01_Reg(Nreg,Ntime),psw02_Reg(Nreg,Ntime),psw03_Reg(Nreg,Ntime))

      Allocate (TPp(Nlon,Nlat,Ntime),TPp_Reg(Nreg,Ntime))

      psw1=Miss; psw2=Miss; psw3=Miss; psw01=Miss; psw02=Miss; psw03=Miss; 
      psw1_Reg=Miss; psw2_Reg=Miss; psw3_Reg=Miss; 
      psw01_Reg=Miss; psw02_Reg=Miss; psw03_Reg=Miss; 
      TPp=Miss; TPp_Reg=Miss

      print*,"In Calculating_SM_Percentile..."
      print*,"Nlon,Nlat,Ntime,Nreg=",Nlon,Nlat,Ntime,Nreg
      Do ilon=1,Nlon
         Do ilat=1,Nlat
            if (Regs(ilon,ilat)>0) then
print*,ilon,ilat
               X(1:Ntime)=sw1(ilon,ilat,1:Ntime)
print*,X(1:Ntime);pause
               Call Sw_percentile (X,Y,Nwd); psw1(ilon,ilat,1:Ntime)=Y(1:Ntime)
pause
               X(1:Ntime)=sw2(ilon,ilat,1:Ntime)
               Call Sw_percentile (X,Y,Nwd); psw2(ilon,ilat,1:Ntime)=Y(1:Ntime)
pause
               X(1:Ntime)=sw3(ilon,ilat,1:Ntime)
               Call Sw_percentile (X,Y,Nwd); psw3(ilon,ilat,1:Ntime)=Y(1:Ntime)
pause
               X(1:Ntime)=sw01(ilon,ilat,1:Ntime)
               Call Sw_percentile (X,Y,Nwd); psw01(ilon,ilat,1:Ntime)=Y(1:Ntime)
               X(1:Ntime)=sw02(ilon,ilat,1:Ntime)
               Call Sw_percentile (X,Y,Nwd); psw02(ilon,ilat,1:Ntime)=Y(1:Ntime)
               X(1:Ntime)=sw03(ilon,ilat,1:Ntime)
               Call Sw_percentile (X,Y,Nwd); psw03(ilon,ilat,1:Ntime)=Y(1:Ntime)

               X(1:Ntime)=TP(ilon,ilat,1:Ntime)    ! Total Precipitation
               Call Sw_percentile (X,Y,Nwd); TPp(ilon,ilat,1:Ntime)=Y(1:Ntime)
            endif
         Enddo
      Enddo

      Do ireg=1,Nreg
         X(1:Ntime)=sw1_Reg(ireg,1:Ntime)
         Call Sw_percentile (X,Y,Nwd); psw1_Reg(ireg,1:Ntime)=Y(1:Ntime)
         X(1:Ntime)=sw2_Reg(ireg,1:Ntime)
         Call Sw_percentile (X,Y,Nwd); psw2_Reg(ireg,1:Ntime)=Y(1:Ntime)
         X(1:Ntime)=sw3_Reg(ireg,1:Ntime)
         Call Sw_percentile (X,Y,Nwd); psw3_Reg(ireg,1:Ntime)=Y(1:Ntime)

         X(1:Ntime)=sw01_Reg(ireg,1:Ntime)
         Call Sw_percentile (X,Y,Nwd); psw01_Reg(ireg,1:Ntime)=Y(1:Ntime)
         X(1:Ntime)=sw02_Reg(ireg,1:Ntime)
         Call Sw_percentile (X,Y,Nwd); psw02_Reg(ireg,1:Ntime)=Y(1:Ntime)
         X(1:Ntime)=sw03_Reg(ireg,1:Ntime)
         Call Sw_percentile (X,Y,Nwd); psw03_Reg(ireg,1:Ntime)=Y(1:Ntime)

         X(1:Ntime)=TP_Reg(ireg,1:Ntime)    ! Total Precipitation
         Call Sw_percentile (X,Y,Nwd); TPp_Reg(ireg,1:Ntime)=Y(1:Ntime)
      Enddo
   End Subroutine Calculating_SM_Percentile
   !----------------------------------------------------------------------------------------------
   Subroutine Sw_percentile (X,Y,Nwd)
      ! Convert soil moisture data from m3/m3 to percentile (0-100)
      ! X: Chuoi SL nhieu nam
      ! Y: Gia tri phan vi cua X
      Integer :: N,i,j,k, nwd, nwd2, N1, Nyy, N60y
      Integer :: Nt(N365)
      Real :: X(Ntime), Y(Ntime), x_value !, X1(N365,N365*Nwd), XX(N365*Nwd), x_value
      Real, Allocatable :: X1(:,:), XX(:)

      Nyy=Ntime/N365   ! So Nam co 365 ngay
      N60y=Nyy*N365    ! So ngay cua tat ca cac nam co du 365 ngay

      Nwd2=Nwd/2; N1=Ntime/N365*Nwd     ! Nyear=Ntime/N365
      Allocate (X1(N365,N1), XX(N1))
      ! Tao mang SL do dai N1=Nyear*Nwd [So nam (Nyear) x So phan tu trong cua so (Nwd)] cua tung ngay trong nam (N365 ngay)
      Do i=1,N365       ! Trich SL cua so nwd cho ngay thu i trong nam
         k=0
         !Do it=i,Ntime,N365         ! it ung voi ngay thu i tung nam
         !   if (it>N60y) exit       ! Chi tinh trong giai doan so nam du 365 ngay
         Do it=i,N60y,N365
            Do j=it-nwd2,it+nwd2    ! Lay Nwd ngay ca truoc va sau ngay thu it
               if (j>0.AND.j<N60y) then ! Ntime) then
                  k=k+1; X1(i,k)=X(j)
!                  print*,"i,k,it,j,X1(i,k)=",i,k,it,j,X1(i,k)
               endif
            Enddo
         Enddo
         Nt(i)=k  ! Nt(i)<=Nyear*Nwd
!print*,"Day ",i," Nt(",i,")=",Nt(i)
      Enddo
      ! 
      Do i=1,Ntime       ! Convert gia tri cua mang X sang gia tri phan vi (0-100) khi so sanh chung voi chuoi X1 o tren
         x_value=X(i)  
         if (x_value==Miss) then
            Y(i)=Miss
         else
            j=mod(i,N365); if (j==0) j=N365; N=Nt(j)   ! Ngay i ung voi ngay j trong nam. N la do dai chuoi X1 cua ngay thu j (1-365) trong nam
            Do k=1,N; XX(k)=X1(j,k); Enddo             ! Chuoi X1 cua ngay thu j trong nam
            Y(i)=Percentile_rank (XX,N,x_value, Miss)  ! Y(i) la gia tri phan vi cua X(i) trong X1
         endif
print*,"==>",i,x_value,Y(i)
      Enddo
   End Subroutine Sw_percentile
   !----------------------------------------------------------------------------------------------
   Subroutine Read_Regs (Fname)
      Character (Len=*) :: Fname
      Integer :: ilon,ilat,Nlo,Nla, i1, i2,i3,i4,i5,i6,i7,i8,i9

      Call rd_ncdf_2d (Fname,Lon,Lat,Nlon,Nlat,Regs,"reg")
   End Subroutine Read_Regs
   !----------------------------------------------------------------------------------------------
   Subroutine Read_sw_data
      Integer :: Nt, it, j, iday, ilon,ilat, ifile, Nobs
      Integer :: Nlo, Nla
      Real, Allocatable :: x2d(:,:), xlon(:),xlat(:),xtime(:),Y1d(:)
      Real :: tmp

      Fname="/home/tanpv/Data/ERA5/sw_19610101-20230612_daily.nc" !era5_1961-2020_daily.nc"              ! Daily data
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,sw1,"sw1")
print*,nlo,nla,nt;pause
print*,sw1(:,:,nt-100:nt);pause
      Deallocate (xlon,xlat,xtime)
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,sw2,"sw2")
      Deallocate (xlon,xlat,xtime)
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,sw3,"sw3")
      Deallocate (xlon,xlat,xtime)

      Fname="/home/tanpv/Data/ERA5/sw_19610101-20230612_pentad_daily.nc" !era5_1961-2020_pentad_daily.nc"       ! Daily pentad data (trung binh truot 5 ngay)
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,sw01,"sw1")
      Deallocate (xlon,xlat,xtime)
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,sw02,"sw2")
      Deallocate (xlon,xlat,xtime)
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,sw03,"sw3")
      Deallocate (xlon,xlat,xtime)

      Fname="/home/tanpv/Data/ERA5/sw_19610101-20230612_daily_areal_mean.nc" !era5_1961-2020_daily_areal_mean.nc"   ! Daily data
      Call rd_ncdf_2d (Fname,xlon,xtime,nlo,nt,sw1_Reg,"sw1")
      Deallocate (xlon,xtime)
      Call rd_ncdf_2d (Fname,xlon,xtime,nlo,nt,sw2_Reg,"sw2")
      Deallocate (xlon,xtime)
      Call rd_ncdf_2d (Fname,xlon,xtime,nlo,nt,sw3_Reg,"sw3")
      Deallocate (xlon,xtime)

      Fname="/home/tanpv/Data/ERA5/sw_19610101-20230612_pentad_daily_areal_mean.nc" !era5_1961-2020_pentad_daily_areal_mean.nc"   ! Daily pentad data (trung binh truot 5 ngay)
      Call rd_ncdf_2d (Fname,xlon,xtime,nlo,nt,sw01_Reg,"sw1")
      Deallocate (xlon,xtime)
      Call rd_ncdf_2d (Fname,xlon,xtime,nlo,nt,sw02_Reg,"sw2")
      Deallocate (xlon,xtime)
      Call rd_ncdf_2d (Fname,xlon,xtime,nlo,nt,sw03_Reg,"sw3")
      Deallocate (xlon,xtime)

      Ntime=Nt

   End Subroutine Read_sw_data
   !----------------------------------------------------------------------------------------------
   Subroutine Write_data_to_file (Fname,Fnm1) ! ,X4d,X1d,Nlon,Nlat,Nt,Nfcs)
      Character (Len=*) :: Fname,Fnm1         ! Daily and Moving Pentad
      ! Dat(Nlon,Nlat,Nfiles,Ndays)
      Integer :: it,ifile,ireg,Ndays_60y
      Real, Allocatable :: X3d(:,:,:),xday(:),X2d(:,:),xreg(:)
      Character (Len=50) :: t_unit, lon_nm, lat_nm, z_nm, t_nm
      Character (Len=10) :: var_nm(Nfiles)=(/"pr","t2m","tx","tm","eva","pet","sw1","sw2","sw3"/)
      Character (Len=10) :: var_unit
      Character (Len=150) :: var_ln(8)=(/"Daily volumetric soil water percemtile layer 1 (0 - 7cm) [%]", &
                                         "Daily volumetric soil water percemtile layer 2 (7 - 28cm) [%]", &
                                         "Daily volumetric soil water percemtile layer 3 (28 - 100cm) [%]", &
                                         "Pentad Daily Total Precipitation percemtile [%]", &
                                         "Pentad Daily 2m-Temperature [C]", &
                                         "Pentad Daily Maximum Temperature [C]", &
                                         "Pentad Daily Minimum Temperature [C]", &
                                         "Pentad Daily Evaporation [mm/day]"/)

      Allocate (xday(Ntime),xreg(Nreg))
      Do ireg=1,Nreg; xreg(ireg)=ireg; Enddo
      Do it=1,Ntime;  xday(it)=it-1; Enddo

      var_unit="%"; t_unit="days since 1961-01-01 ";  lon_nm="lon"; lat_nm="lat"; t_nm="day"
      ! So lieu SW ngay da convert sang phan vi (percentile) 

      print*,"Writing sw01..."
      Call wr_ncdf_3d (Fname,Lon,Lat,xday,Nlon,Nlat,Ntime,psw01,"sw01",t_unit,"Pentad "//trim(var_ln(1)),var_unit,lon_nm,lat_nm,t_nm,Miss)
      print*,"Writing sw02..."
      Call wr_ncdf_3d_apd (Fname,psw02,"sw02",Nlon,Nlat,Ntime,"Pentad "//trim(var_ln(2)),var_unit,lon_nm,lat_nm,t_nm,Miss)
      print*,"Writing sw03..."
      Call wr_ncdf_3d_apd (Fname,psw03,"sw03",Nlon,Nlat,Ntime,"Pentad "//trim(var_ln(3)),var_unit,lon_nm,lat_nm,t_nm,Miss)

      ! Areal mean of daily data:
      Do ireg=1,Nreg; xreg(ireg)=ireg; Enddo

      Call wr_ncdf_2d (Fnm1,xreg,xday,Nreg,Ntime,psw01_Reg,"sw01","Areal mean of pentad "//trim(var_ln(1)),var_unit,"reg","day",Miss)
      Call wr_ncdf_2d_apd (Fnm1,psw02_Reg,"sw02",Nreg,Ntime,"Areal mean of pentad "//trim(var_ln(2)),var_unit,"reg","day",Miss)
      Call wr_ncdf_2d_apd (Fnm1,psw03_Reg,"sw03",Nreg,Ntime,"Areal mean of pentad "//trim(var_ln(3)),var_unit,"reg","day",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/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