Use NetCDF

   Real, Parameter :: Miss=-99.0
   Integer, Parameter :: IMiss=-99, Nreg=7, Pentad=5
!   Integer, Parameter :: Period=1, N5y=365   ! 365 days x 1 years (period)
   Integer, Parameter :: Period=5, N5y=1825  ! 1825=365 days x 5 years (period)
   Integer, Parameter :: Ncha=5   ! b(1), number of FD; b(2), mean duration of FD; b(3), mean severity of FD; b(4), mean speed of FD; b(5), FD intensification rate
   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   ! (Nlon,Nlat,Ntime) 
   Real, Allocatable, Dimension (:,:)   :: sw1_Reg,sw2_Reg,sw3_Reg ! (Nreg,Ntime)
   
   Real, Allocatable, Dimension (:,:,:,:) :: psw5      ! (Nlon,Nlat,Nyrs,N5y) - Tinh cho tung thoi ky 5 nam
   Real, Allocatable, Dimension (:,:,:,:) :: FDCha5  ! (Nlon,Nlat,Ncha,Nyrs)
   Real, Allocatable, Dimension (:,:,:)   :: psw5_Reg  ! (Nreg,Nyrs,N5y)
   Real, Allocatable, Dimension (:,:,:)   :: FDCha5_Reg ! (Nreg,Ncha,Nyrs)
   

   Real, Allocatable, Dimension (:,:,:) :: psw1   ! (Nlon,Nlat,Ntime) 
   Real, Allocatable, Dimension (:,:,:) :: FDevt  ! drt(Nlon,Nlat,Ntime) - 1-drought onset, 2-drought recovery, 0-no drought
   Real, Allocatable, Dimension (:,:,:) :: FDCha  ! (Nlon,Nlat, Ncha)
   Real, Allocatable, Dimension (:,:)   :: psw1_Reg ! (Nreg,Ntime)
   Real, Allocatable, Dimension (:,:)   :: FDevt_Reg  ! drt(Nreg,Ntime)
   Real, Allocatable, Dimension (:,:)   :: FDCha_Reg  ! (Nreg, Ncha)

   Real, Allocatable, Dimension (:,:,:,:) :: X4d
   Integer, Allocatable, Dimension (:) :: date_m1   ! (Nfcs)
   Integer :: Nlon, Nlat, Nlev, Ntime, Nyrs, Ndays, Npentad, Nlon0, Nlat0    ! Nfcs = So lan du bao, date_fcs = yyyymmdd cua ngay DB

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

   OutDir="/home/tanpv/Progs/Flash_Dr/NC/"
   Call Read_Regs ("/home/tanpv/Data/ERA5/Mask_VN_0.25.nc") !VN.grid_0.25.nc")

   Call Read_sw_percentile_data
   Call Calculating_FD_Chars
   if (Period==5) then
      Call Write_data_to_file (trim(OutDir)//"FD_Char_1961-2020_5yrs_Daily.nc",trim(OutDir)//"FD_Char_1961-2020_5yrs_Daily_Reg.nc")
   else
      Call Write_data_to_file (trim(OutDir)//"FD_Char_1961-2020_1yrs_Daily.nc",trim(OutDir)//"FD_Char_1961-2020_1yrs_Daily_Reg.nc")
   endif
   CONTAINS
   !----------------------------------------------------------------------------------------------
   !Include "FL_Dr_1.f90"
   Include "Detect_event.f90"
   Include "Sub_Spell.f90"
   !----------------------------------------------------------------------------------------------
   Subroutine FD_Chars (X,N,b)
      Real, Parameter :: Th40=40., Th20=20., Speed=5.
      Integer, Parameter :: Lmin=15   ! Do dai toi thieu Lmin days
      Integer :: N, it, i1, i2, i, j, k, k1, cnt, ic, kid
      Real :: X(N), Y40(N), Y20(N), b(Ncha)
      Real :: M40(N), IMAX40, Spells40(N), N_Spell40, Start40(N), End40(N)
      Real :: M20(N), IMAX20, Spells20(N), N_Spell20, Start20(N), End20(N)
      Integer :: N40,N20,I40,I20, ist, ien
      Real :: Score, Spd, Spdx, Spd1, swmin, osev(n), odur(n), ospd(n)
      Real :: RDi, RD(n) ! FD intensification rate (RD) = Sum[ (sw(i-1)-sw(i)) ]/Duration
      ! It's modified by: RD = [(sw(ist-1)-swmin)/sw(ist-1)*100]/[it(sw=swmin)-ist+1]
      ! (Mukherjee, S., & Mishra, A. K. (2022). A multivariate flash drought indicator for identifying global hotspots and associated climate controls. 
      !                                         Geophysical Research Letters, 49, e2021GL096804. https://doi. org/10.1029/2021GL096804
      Logical :: Less,flag

      Less=.True.; Score=1  ! sw <=40
      Call Events (X,Y40,N,Th40,Score,Less,Miss)
      Score=1   ! sw<=20
      Call Events (X,Y20,N,Th20,Score,Less,Miss)

      Call NSpell_and_Leng (Y40,N,M40,Imax40,Spells40,N_Spell40,Start40,End40,Miss)
      Call NSpell_and_Leng (Y20,N,M20,Imax20,Spells20,N_Spell20,Start20,End20,Miss)
      N40=N_Spell40; N20=N_Spell20
      !Do i=1,N20
      !   if (Int(Spells20(i))<Lmin) then    ! Bo qua cac dot X<20 it hon Lmin ngay
      !      i1=Int(Start20(i)); i2=Int(End20(i))
      !      Y20(i1:i2)=0
      !   endif
      !Enddo
      Spdx=0; cnt=0     ! So FD
!print*,"Begin FD_Chars ....................................."
      Loop1: Do i=1,N40
         if (Int(Spells40(i))<Lmin) Cycle Loop1
         ist=Int(Start40(i)); ien=Int(End40(i))
!print"(50F5.1)",X(ist-3:ien+3)
         j=ist; swmin=100.; flag=.true.
         Loop2: Do
            j=j+1; if (j>ien) exit
            if (Y20(j)/=1) Cycle Loop2
            if (j>1.AND.X(j)/=Miss.AND.X(j)<X(j-1).AND.X(j)<swmin) then; swmin=X(j); kid=j; endif
            
!            print*,X(j-1),X(j)
            !Spd=(X(ist)-X(j))/Real(j-ist)*Pentad   ! Difference of sw at initial and 20% point
            Spd=(Th40-X(j))/Real(j-ist)*Pentad
            if (Spdx<Spd) Spdx=Spd
            if (Spd>=Speed.AND.(j-ist+1)>Pentad) then
               if (flag) then; cnt=cnt+1; flag=.false.; endif
               !Spdx=Spd;
               k=j
               Loop4: Do
                  k=k+1; if (Y20(k)/=1) then; k=k-1; exit; endif   ! Until X>20%
                  if (k<=ien.AND.X(k)/=Miss.AND.X(k)<X(k-1).AND.X(k)<swmin) then; swmin=X(k); kid=k; endif
               Enddo Loop4

               k1=ien+1
               Loop3: Do
                  k1=k1-1; if (k1==k+1) exit
                  if (Y20(k1)==1) exit
                  if (Y20(k1)==0) Cycle Loop3
               Enddo Loop3
               Do it=k,k1
               !!   if (X(it)/=Miss.AND.X(it)<swmin) then; swmin=X(it); RDi=it; endif
                  
!                  print*,swmin
                  !Spd1=(X(ist)-X(it))/Real(it-ist)*Pentad
                  Spd1=(Th40-X(it))/Real(it-ist)*Pentad
                  if (Spd1>Spdx) Spdx=Spd1
               Enddo
               k=k1

               !Loop3: Do
               !   k=k+1; if (Y20(k)==1) Cycle Loop3; 
                  !if (Y20(k)/=1) exit
               !   if (Y20(k)/=1) then
               !      k1=k; ic=0
               !      Loop4: Do
               !         k1=k1+1; if (k1>=ien) exit
               !         if (Y20(k1)/=1) then; ic=ic+1; Cycle Loop4; endif
               !      Enddo Loop4
               !   endif
               !   if (ien-k1>10) k=k1; 
               !   if (Y20(k)/=1) exit
               !   if (k>ien) exit
               !   if (X(k)/=Miss.AND.X(k)<swmin) swmin=X(k)
               !   Spd1=(X(ist)-X(k))/Real(k-ist)*Pentad
               !   if (Spd1>Spdx) Spdx=Spd1
               !Enddo Loop3
               odur(cnt)=k-ist+1     ! Duration in days
               osev(cnt)=0
               Do it=ist,k
                  osev(cnt)=osev(cnt)+(Th40-X(it))
               Enddo
               ospd(cnt)=Spdx
               RDi=0
               Do it=ist,kid
                  if (it==ist) then
                     RDi=RDi+(Th40-X(it))/Th40*100.
                  else
                     RDi=RDi+(X(it-1)-X(it))/X(it-1)*100.
                  endif
               Enddo
               if (kid>ist) then
                  RD(cnt)=RDi/Real(kid-ist+1)
               else
                  RD(cnt)=Miss
               endif
               !!if (ist>1) then 
               !!   RD(cnt)=(X(ist-1)-swmin)/X(ist-1)*100.0/(RDi-real(ist)+1.)
               !!else
               !!   RD(cnt)=(X(ist)-swmin)/X(ist)*100.0/(RDi-real(ist))
               !!endif
!print*,RDi
            endif
         Enddo Loop2
         if (odur(cnt)<Lmin) then; if (cnt>1) then; cnt=cnt-1; else; cnt=0; endif; endif
!print*,"Spell: ",i,ist,ien,k
!print"(50F5.1)",X(ist-3:ien+3)
!print"(50F3.0)",Y40(ist-3:ien+3)
!print"(50F3.0)",Y20(ist-3:ien+3)
!print*,cnt,odur(cnt),osev(cnt),ospd(cnt)

!pause
      Enddo Loop1
      if (cnt==0) then 
         b=0.
      else
         ! Trung binh nam:
         b(1)=cnt/Real(Period);                            ! Number of FD
         b(2)=Xsum(odur,cnt,Miss)/real(cnt)/Real(Period);  ! Mean duration
         b(3)=Xsum(osev,cnt,Miss)/real(cnt)/Real(Period);  ! Mean severity
         b(4)=Xsum(ospd,cnt,Miss)/real(cnt)/Real(Period);  ! Mean of maximum decline speed
         b(5)=Xsum(RD,cnt,Miss)/real(cnt)/Real(Period)     ! Mean intensification rate
      endif
!print*,"===> ",b; pause
   End Subroutine FD_Chars
   !----------------------------------------------------------------------------------------------
   Subroutine Calculating_FD_Chars
      Integer :: ilon, ilat, it, ireg
      Real :: X(Ntime), drt(Ntime), b(Ncha)
      ! b(1), number of FD; b(2), mean duration of FD; b(3), mean severity of FD; b(4), mean speed of FD

      Allocate (FDCha5(Nlon,Nlat,Ncha,Nyrs))
      Allocate (FDCha5_Reg(Nreg,Ncha,Nyrs))

      FDCha5=Miss; FDCha5_Reg=Miss

      Do ilon=1,Nlon
         Do ilat=1,Nlat
            if (Regs(ilon,ilat)>0) then
               Do it=1,Nyrs
                  X(1:Ntime)=psw5(ilon,ilat,it,1:Ntime)
                  !Call flashdrought(Ntime,X,b,drt)
                  Call FD_Chars (X,Ntime,b)
                  !b(1) = b(1)/real(Pentad)                ! number of FD per year
                  FDCha5(ilon,ilat,1:Ncha,it) =b(1:Ncha)
               Enddo
            endif
         Enddo
      Enddo

      Do ireg=1,Nreg
         Do it=1,Nyrs
            X(1:Ntime)=psw5_Reg(ireg,it,1:Ntime)
            !Call flashdrought(Ntime,X,b,drt)
            Call FD_Chars (X,Ntime,b)
            !b(1) = b(1)/real(Pentad)                ! number of FD per year
            FDCha5_Reg(ireg,1:Ncha,it) =b(1:Ncha)
         Enddo
      Enddo
   End Subroutine Calculating_FD_Chars
   !----------------------------------------------------------------------------------------------
   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_percentile_data
      Integer :: Nt, it, j, iday, ilon,ilat, ifile, Nobs, i1,i2
      Integer :: Nlo, Nla
      Real, Allocatable :: xlon(:),xlat(:),xtime(:),Y1d(:)
      Real :: tmp

      Fname=trim(OutDir)//"percentile_sw_1961-2020_daily.nc"   !  gridpoint
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,sw1,"sw01")
      Deallocate (xlon,xlat,xtime)
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,sw2,"sw02")
      Deallocate (xlon,xlat,xtime)
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,sw3,"sw03")
      Deallocate (xlon,xlat,xtime)

      ! 3-layer mean:
      Allocate (psw1(nlo,nla,nt)); psw1=Miss
      Do ilon=1,Nlon
         Do ilat=1,Nlat
            Do ireg=1,Nreg
               if (Regs(ilon,ilat)>0) then
                  Do it=1,Nt
                     j=0; tmp=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) psw1(ilon,ilat,it)=tmp/real(j)
                  Enddo
               endif
            Enddo
         Enddo
      Enddo
!      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,psw1,"swp")
!      Deallocate (xlon,xlat,xtime)

      Fname=trim(OutDir)//"percentile_sw_1961-2020_daily_Reg.nc"  ! Regional mean
      Call rd_ncdf_2d (Fname,xlon,xtime,nlo,nt,sw1_Reg,"sw01")
      Deallocate (xlon,xtime)
      Call rd_ncdf_2d (Fname,xlon,xtime,nlo,nt,sw2_Reg,"sw02")
      Deallocate (xlon,xtime)
      Call rd_ncdf_2d (Fname,xlon,xtime,nlo,nt,sw3_Reg,"sw03")
      Deallocate (xlon,xtime)

      ! 3-layer mean:
      Allocate (psw1_Reg(Nreg,Nt)); psw1_Reg=Miss
      Do ireg=1,Nreg
         Do it=1,Nt
            j=0; tmp=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) psw1_Reg(ireg,it)=tmp/real(j)
         Enddo
      Enddo
!      Call rd_ncdf_2d (Fname,xlon,xtime,nlo,nt,psw1_Reg,"swp")
!      Deallocate (xlon,xtime)

      Nyrs=Nt/N5y
      Ntime=N5y

      Allocate (psw5(Nlon,Nlat,Nyrs,N5y),psw5_Reg(Nreg,Nyrs,N5y))
      Do ilon=1,Nlon
         Do ilat=1,Nlat
            Do j=1,Nyrs
               i1=(j-1)*N5y+1; i2=j*N5y
               psw5(ilon,ilat,j,1:N5y)=psw1(ilon,ilat,i1:i2)
            Enddo
         Enddo
      Enddo
      Do ireg=1,Nreg
         Do j=1,Nyrs
            i1=(j-1)*N5y+1; i2=j*N5y
            psw5_Reg(ireg,j,1:N5y)=psw1_Reg(ireg,i1:i2)
         Enddo
      Enddo
   End Subroutine Read_sw_percentile_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,Ntime
      Real, Allocatable :: X3d(:,:,:),xday(:),X2d(:,:),xreg(:),xcha(:)
      Character (Len=50) :: t_unit, lon_nm, lat_nm, z_nm, t_nm
      Character (Len=10) :: var_unit
      Character (Len=150) :: var_ln(3)=(/"FD chars: b(1), number of FD per year; b(2), mean duration of FD; b(3), mean severity of FD;"// &
                                               " b(4), mean speed (%/pentad) or intensification rate (%) of FD", &
                                         "Daily volumetric soil water percemtile layer 2 (7 - 28cm)", &
                                         "Daily volumetric soil water percemtile layer 3 (28 - 100cm)"/)

!      Allocate (FDCha5(Nlon,Nlat,Ncha,Nyrs))
!      Allocate (FDCha5_Reg(Nreg,Ncha,Nyrs))
      Ntime=Nyrs
      Allocate (xday(Ntime),xreg(Nreg),xcha(Ncha))
      Do ireg=1,Nreg; xreg(ireg)=ireg; Enddo
      Do it=1,Ntime;  xday(it)=(it-1)*Period; Enddo
      Do it=1,Ncha;   xcha(it)=it; Enddo

      var_unit="-"; t_unit="years since 1961-01-01 ";  lon_nm="lon"; lat_nm="lat"; z_nm="cha"; t_nm="year"
      Call wr_ncdf_4d (Fname,Lon,Lat,xcha,xday,Nlon,Nlat,Ncha,Nyrs,FDCha5,"fd", &
                       t_unit, trim(var_ln(1)), var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss)
!      Call wr_ncdf_3d (Fname,Lon,Lat,xcha,Nlon,Nlat,Ncha,FDCha,"fd",t_unit,trim(var_ln(1)),var_unit,lon_nm,lat_nm,t_nm,Miss)
!      Call wr_ncdf_3d_apd (Fname,psw2,"sw2",Nlon,Nlat,Ntime,trim(var_ln(2)),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)

      ! Areal mean of daily data:
      Do ireg=1,Nreg; xreg(ireg)=ireg; Enddo
      lon_nm="reg"; lat_nm="cha"; t_nm="year"
      Call wr_ncdf_3d (Fnm1,xreg,xcha,xday,Nreg,Ncha,Nyrs,FDCha5_Reg,"fd",t_unit,trim(var_ln(1)),var_unit,lon_nm,lat_nm,t_nm,Miss)
!      Call wr_ncdf_2d (Fnm1,xreg,xcha,Nreg,Ncha,FDCha_Reg,"fd",trim(var_ln(1)),var_unit,"reg","cha",Miss)
!      Call wr_ncdf_2d_apd (Fnm1,psw2_Reg,"sw2",Nreg,Ntime,"Areal mean of "//trim(var_ln(2)),var_unit,"reg","day",Miss)
!      Call wr_ncdf_2d_apd (Fnm1,psw3_Reg,"sw3",Nreg,Ntime,"Areal mean of "//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