Use NetCDF

   Real, Parameter :: Miss=-99.0
   Integer, Parameter :: IMiss=-99, Nfiles=9, Nreg=7, Pentad=5, N365=73, Nwd=3  ! Nwd=Window ata
   Character (Len=25), Parameter :: Fnames(Nfiles)=(/"VN_tp.nc","VN_t2m.nc","VN_tx.nc","VN_tm.nc","VN_eva.nc","VN_pet.nc","VN_sw1.nc","VN_sw2.nc","VN_sw3.nc"/)
   Character (Len=10), Parameter :: vname(Nfiles)=(/"tp","t2m","mx2t","mn2t","e","pev","swvl1","swvl2","swvl3"/)
   !Character (Len=25), Parameter :: Fnames(Nfiles)=(/"tp_daily.nc","t2m_daily.nc","tx_daily.nc","tm_daily.nc","eva_daily.nc", &
   !                                                  "pet_daily.nc","sw1_daily.nc","sw2_daily.nc","sw3_daily.nc"/)
   !Character (Len=10), Parameter :: vname(Nfiles)=(/"var228","t2m","mx2t","mn2t","var182","var251","swvl1","swvl2","swvl3"/)
   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 (:,:,:) :: sw13p                        ! (Nlon,Nlat,Npentad) - Trung binh pentad cua sw1-sw3
   Real, Allocatable, Dimension (:,:)   :: sw1_Reg,sw2_Reg,sw3_Reg,sw01_Reg,sw02_Reg,sw03_Reg ! (Nreg,Ntime)
   Real, Allocatable, Dimension (:,:)   :: sw13p_Reg                    ! (Nreg,Npentad) - Trung binh pentad cua sw1-sw3 lay trung binh vung
   
   Real, Allocatable, Dimension (:,:,:) :: psw1,psw2,psw3,psw01,psw02,psw03   ! (Nlon,Nlat,Ntime) 
   Real, Allocatable, Dimension (:,:,:) :: psw13p                       ! (Nlon,Nlat,Npentad) - percentile cua Trung binh pentad cua sw1-sw3

   Real, Allocatable, Dimension (:,:)   :: psw1_Reg,psw2_Reg,psw3_Reg,psw01_Reg,psw02_Reg,psw03_Reg ! (Nreg,Ntime)
   Real, Allocatable, Dimension (:,:)   :: psw13p_Reg                   ! (Nreg,Npentad) - percentile cua psw13p

   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
!   Real :: XX(30),YY(30)
!   Data XX /22,69,49,78,41,39,24,48,35,79,31,72,94,37,90,10,25,14,92,24,24,41,54,55,60,95,75,42,17,75/

!   Ntime=30
   
!   print"(50F4.0)",XX(:)
!   print"(50F4.0)",YY(:)
!   Call Sw_percentile (XX,YY,5)
!pause

!      Fname="tp.nc"
!         Call GetTime (Fname, Nt, factor,offset,Missing,fill_value)
!         Call rd_ncdf_4d (Fname,lon,lat,lev,xtime,Nlon,Nlat,Nlev,Ntime,X4d,"tp")
!      Where (X4d==Missing) X4d=Miss
!      Where (X4d/=Miss) X4d=X4d * factor + offset ! unpacked_value = packed_value * scale_factor + add_offset

!print*,Nlon,Nlat,Nlev,Ntime

!print*,X4d(:,:,2,:); pause
!print*,X4d(:,:,2,:)

   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)//"Pentad_sw_1961-2020.nc",trim(OutDir)//"Pentad_sw_1961-2020_Reg.nc")
   CONTAINS
   !----------------------------------------------------------------------------------------------
   Subroutine Calculating_SM_Percentile
      Integer :: ilon, ilat, it, ireg, Ntime
      Real, Allocatable :: X(:), Y(:)

      Ntime=Npentad
print*,"Ntime,Nwd=",Ntime,Nwd;pause
      Allocate (psw1(Nlon,Nlat,Ntime),X(Ntime),Y(Ntime))
      Allocate (psw1_Reg(Nreg,Ntime))

      psw1=Miss
      psw1_Reg=Miss

      Do ilon=1,Nlon
         Do ilat=1,Nlat
            if (Regs(ilon,ilat)>0) then
print*,ilon,ilat
               X(1:Ntime)=sw13p(ilon,ilat,1:Ntime)
               Call Sw_percentile (X,Y,Ntime,N365,Nwd); psw1(ilon,ilat,1:Ntime)=Y(1:Ntime)
            endif
         Enddo
      Enddo

      Do ireg=1,Nreg
         X(1:Ntime)=sw13p_Reg(ireg,1:Ntime)
         Call Sw_percentile (X,Y,Ntime,N365,Nwd); psw1_Reg(ireg,1:Ntime)=Y(1:Ntime)
      Enddo
   End Subroutine Calculating_SM_Percentile
   !----------------------------------------------------------------------------------------------
   Subroutine Sw_percentile (X,Y,Ntime,N365,Nwd)
      ! Convert soil moisture data from m3/m3 to percentile (0-100)
      Integer :: N,i,j,k, nwd, nwd2, N1
      Integer :: Nt(N365)
      Real :: X(Ntime), Y(Ntime), x_value !, X1(N365,N365*Nwd), XX(N365*Nwd), x_value
      Real, Allocatable :: X1(:,:), XX(:)

      Nwd2=Nwd/2; N1=Ntime/N365*Nwd
      Allocate (X1(N365,N1), XX(N1))
      Do i=1,N365       ! Trich SL cua so nwd
         k=0
         Do it=i,Ntime,N365
            Do j=it-nwd2,it+nwd2
               if (j>0.AND.j<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)<=Ntime*Nwd
      Enddo
      ! 
      Do i=1,Ntime
         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)
            Do k=1,N; XX(k)=X1(j,k); Enddo
            Y(i)=Percentile_rank (XX,N,x_value, Miss)
         endif
      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(:,:), X3d(:,:,:), xlon(:),xlat(:),xtime(:),Y1d(:)
      Real :: X1d(Pentad), tmp

      Fname="/home/tanpv/Data/ERA5/era5_1961-2020_daily.nc"
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,sw1,"sw1")
      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/era5_1961-2020_daily_areal_mean.nc"
      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)

      Ntime=Nt; Npentad=Ntime/Pentad
      
      Allocate (sw13p(Nlon,Nlat,Npentad),sw13p_Reg(Nreg,Npentad))
      Allocate (X3d(Nlon,Nlat,Ntime),X2d(Nreg,Ntime))
      sw13p=Miss; sw13p_Reg=Miss

      Do ilon=1,Nlon
         Do ilat=1,Nlat
            Do ireg=1,Nreg
               if (Regs(ilon,ilat)>0) then
                  Do it=1,Ntime
                     tmp=0.0; j=0
                     if (sw1(ilon,ilat,it)/=Miss) then; j=j+1; tmp=tmp+sw1(ilon,ilat,it); endif
                     if (sw2(ilon,ilat,it)/=Miss) then; j=j+1; tmp=tmp+sw2(ilon,ilat,it); endif
                     if (sw3(ilon,ilat,it)/=Miss) then; j=j+1; tmp=tmp+sw3(ilon,ilat,it); endif
                     if (j>0) then; X3d(ilon,ilat,it)=tmp/real(j); else; X3d(ilon,ilat,it)=Miss; endif
                  Enddo
                  j=0
                  Do it=1,Ntime,Pentad
                     X1d(1:Pentad)=X3d(ilon,ilat,it:it+Pentad-1)
                     j=j+1; sw13p(ilon,ilat,j)=Ave(X1d,Pentad,Miss)
                  Enddo
               endif
            Enddo
         Enddo
      Enddo
      Do ireg=1,Nreg
         Do it=1,Ntime
            tmp=0.0; j=0
            if (sw1_Reg(ireg,it)/=Miss) then; j=j+1; tmp=tmp+sw1_Reg(ireg,it); endif
            if (sw2_Reg(ireg,it)/=Miss) then; j=j+1; tmp=tmp+sw2_Reg(ireg,it); endif
            if (sw3_Reg(ireg,it)/=Miss) then; j=j+1; tmp=tmp+sw3_Reg(ireg,it); endif
            if (j>0) then; X2d(ireg,it)=tmp/real(j); else; X2d(ireg,it)=Miss; endif
         Enddo
         j=0
         Do it=1,Ntime,Pentad
            X1d(1:Pentad)=X2d(ireg,it:it+Pentad-1)
            j=j+1; sw13p_Reg(ireg,j)=Ave(X1d,Pentad,Miss)
         Enddo
      Enddo
   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(4)=(/"Pentad volumetric soil water at Root zone layer (0 - 100cm)", &
                                         "Regional mean of Pentad volumetric soil water at Root zone layer  (0 - 100cm)", &
                                         "Pentad volumetric soil water percemtile at Root zone layer  (0 - 100cm)", &
                                         "Regional mean of Pentad volumetric soil water percemtile at Root zone layer  (0 - 100cm)"/)


      Allocate (xday(Npentad),xreg(Nreg))
      Do ireg=1,Nreg; xreg(ireg)=ireg; Enddo
      Do it=1,Npentad;  xday(it)=(it-1)*Pentad; Enddo

      var_unit="m3/m3"; t_unit="days since 1961-01-01 ";  lon_nm="lon"; lat_nm="lat"; t_nm="pentad"
      Call wr_ncdf_3d (Fname,Lon,Lat,xday,Nlon,Nlat,Npentad,sw13p,"sw",t_unit,trim(var_ln(1)),var_unit,lon_nm,lat_nm,t_nm,Miss)
      var_unit="%";
      Call wr_ncdf_3d_apd (Fname,psw1,"swp",Nlon,Nlat,Npentad,trim(var_ln(3)),var_unit,lon_nm,lat_nm,t_nm,Miss)

!      Call wr_ncdf_3d_apd (Fname,psw3,"sw3",Nlon,Nlat,Ntime,trim(var_ln(3)),var_unit,lon_nm,lat_nm,t_nm,Miss)

!      Call wr_ncdf_3d_apd (Fname,psw01,"sw01",Nlon,Nlat,Ntime,"Pentad "//trim(var_ln(1)),var_unit,lon_nm,lat_nm,t_nm,Miss)
!      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)
!      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
      var_unit="m3/m3"
      Call wr_ncdf_2d (Fnm1,xreg,xday,Nreg,Npentad,sw13p_Reg,"sw",trim(var_ln(2)),var_unit,"reg","pentad",Miss)
      var_unit="%";
      Call wr_ncdf_2d_apd (Fnm1,psw1_Reg,"swp",Nreg,Npentad,trim(var_ln(4)),var_unit,"reg","pentad",Miss)

!      Call wr_ncdf_2d_apd (Fnm1,psw3_Reg,"sw3",Nreg,Ntime,"Areal mean of "//trim(var_ln(3)),var_unit,"reg","day",Miss)

!      Call wr_ncdf_2d_apd (Fnm1,psw01_Reg,"sw01",Nreg,Ntime,"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