! Tinh cac dac trung thong ke cua chi so han
Use netcdf
Implicit None
Integer, Parameter :: Yr1=1986, Yr2=2005, Nmon=12, Nyear=(Yr2-Yr1)+1
Integer, Parameter :: Ntime=Nyear*Nmon, Nreg=7, Nst=131, Nbin=18
Integer, Parameter :: Nlon=40, Nlat=68

Real, Parameter :: Miss = -99.0

Real, Parameter :: Thre_1_PDSI=-2, Thre_2_PDSI=-3, Thre_3_PDSI=-4   ! PDSI: Moderate, Severe, Extreme
Real, Parameter :: BIN_PDSI(Nbin+1)=(/-30,-20,-15,-10,-6,-4,  -3,-2,-1,  0,  1,2,  3,4,6,10,15,20,30/)

Real, Parameter :: Thre_1_SPI=-1, Thre_2_SPI=-1.5, Thre_3_SPI=-2 !  SPI: Moderate, Severe, Extreme
Real, Parameter :: BIN_SPI(Nbin+1) =(/-15.,-10., -6., -4.,-3.,-2.,-1.5,-1.,-0.5,0.,0.5,1.,1.5,2.,3., 4., 6.,10.,15./)
Real :: Thre_1, Thre_2, Thre_3, BIN(Nbin+1)

Real :: Freq(Nreg, Nbin), Freq_Mon(Nreg, Nbin, Nmon), Xbin(Nbin), Xreg(Nreg)    ! Frequency

! Definitions:
Real, Dimension(Nlon,Nlat,Ntime) :: D1, D2, D3    ! (Nlon, Nlat, Ntime); =0/1; D1: Moderate+Severe+Extreme, D2: Severe+Extreme, D3: Extreme
Real, Dimension(Nlon,Nlat,Ntime) :: OD1, OD2, OD3 ! (Nlon, Nlat, Ntime)  =1: No drought event

! Events and Duration, and Time Arrival: (Nlon, Nlat, Ntime), (Nlon, Nlat)
Real, Dimension(Nlon,Nlat,Ntime) :: Dur1, Dur2, Dur3        ! Tan so cac dot
Real, Dimension(Nlon,Nlat)       :: N_Dur1, N_Dur2, N_Dur3  ! Max cua Do dai cac dot

Real, Dimension(Nlon,Nlat,Ntime) :: Arr1, Arr2, Arr3        ! Time Interval Frequency ? 
Real, Dimension(Nlon,Nlat)       :: N_Arr1, N_Arr2, N_Arr3  ! (Nlon, Nlat) Max cua khoang thoi gian giua hai dot

Real, Dimension(Nlon,Nlat,Ntime) :: Spel1, Spel2, Spel3 ! Do dai tung dot cua tung tram
Real, Dimension(Nlon,Nlat,Ntime) :: DS1, DS2, DS3, DI1, DI2, DI3 ! Severity & Intensity of each spell at each grid    <<<<<<<<<
Real, Dimension(Nlon,Nlat)       :: N_Spel1, N_Spel2, N_Spel3 ! So dot trong toan chuoi cua tung tram
Real, Dimension(Nlon,Nlat,Ntime) :: IS1, IS2, IS3   ! Chi so dau cua cac dot trong chuoi
Real, Dimension(Nlon,Nlat,Ntime) :: IE1, IE2, IE3   ! Chi so cuoi cua cac dot trong chuoi

! Regional mean:
Real :: Dur_Reg1(Nreg,Ntime), Arr_Reg1(Nreg,Ntime), NDur_Reg1(Nreg), NArr_Reg1(Nreg)  ! ! Tan so cac dot va Max cua Do dai cac dot cua tung vung
Real :: Dur_Reg2(Nreg,Ntime), Arr_Reg2(Nreg,Ntime), NDur_Reg2(Nreg), NArr_Reg2(Nreg)  ! ! Tan so cac dot va Max cua Do dai cac dot cua tung vung
Real :: Dur_Reg3(Nreg,Ntime), Arr_Reg3(Nreg,Ntime), NDur_Reg3(Nreg), NArr_Reg3(Nreg)  ! ! Tan so cac dot va Max cua Do dai cac dot cua tung vung

Integer :: Max_So_Dot        ! So dot lon nhat cua tat ca cac tram (=Max {Ndot} - Gia tri lon nhat cua N_Spell)
Integer :: Max_Length_Dot    ! Do keo dai dai nhat trong tat ca cac dot cua cac tram (=Max {MaxDot} - Gia tri lon nhat cua Imax)
Integer :: Max_So_Dot_Reg        ! So dot lon nhat cua tat ca cac tram (=Max {Ndot} - Gia tri lon nhat cua N_Spell)
Integer :: Max_Length_Dot_Reg    ! Do keo dai dai nhat trong tat ca cac dot cua cac tram (=Max {MaxDot} - Gia tri lon nhat cua Imax)

Real :: Spel_Reg1(Nreg,Ntime), Spel_Reg2(Nreg,Ntime), Spel_Reg3(Nreg,Ntime)    ! Do dai tung dot cua tung vung
Real, Dimension(Nreg,Ntime) :: DS1_Reg, DS2_Reg, DS3_Reg, DI1_Reg, DI2_Reg, DI3_Reg ! ! Severity & Intensity of each spell for each sub-region    <<<<<<<<<
Real, Dimension(Nreg,Ntime) :: DS1_Reg_Std, DS2_Reg_Std, DS3_Reg_Std, DI1_Reg_Std, DI2_Reg_Std, DI3_Reg_Std   ! Do lech chuan
Real :: NSpel_Reg1(Nreg), NSpel_Reg2(Nreg), NSpel_Reg3(Nreg)                   ! So dot trong toan chuoi cua tung vung

Real :: IS_Reg1(Nreg,Ntime), IS_Reg2(Nreg,Ntime), IS_Reg3(Nreg,Ntime)   ! Chi so dau cua cac dot trong chuoi
Real :: IE_Reg1(Nreg,Ntime), IE_Reg2(Nreg,Ntime), IE_Reg3(Nreg,Ntime)   ! Chi so cuoi cua cac dot trong chuoi

Real :: D1_reg(Nreg,Ntime), D2_Reg(Nreg,Ntime), D3_Reg(Nreg,Ntime)      ! =0/1
Real :: OD1_reg(Nreg,Ntime), OD2_Reg(Nreg,Ntime), OD3_Reg(Nreg,Ntime)   ! No drought events

! Geographic Area Extent:
Real :: Area_Reg1(Nreg,Ntime), Area_Reg2(Nreg,Ntime), Area_Reg3(Nreg,Ntime)  !

Real :: PDSI(Nlon,Nlat,Ntime), PDSI_Reg(Ntime,Nreg) ! Drought Index
Real, Allocatable :: Dr_IDX_Tmp(:,:,:)              ! Luu tam thoi PDSI va SPI khi doc SL
Real :: Lon(Nlon), Lat(Nlat), Reg(Nlon,Nlat)

Real, allocatable, dimension (:)   :: loninfo, latinfo
Real, allocatable, dimension (:,:) :: vietnam ! thong tin cac vung khi hau: B1=2, B2=3,..., N3=8
Integer :: nxlon,nxlat

Character (Len=255) :: Fname, Models
!------------------------------------------------------------------------------
Integer, Parameter :: NFile = 36  ! 6 models + 6 models * 2 Scenarios * 2 Future Periods + 5 Ens + 1 Obs
Character (Len=100) :: Fname_inp(NFile), Fname_Out(NFile)
Namelist /Filenames/ Fname_inp, Fname_Out
Integer :: ifile, iidx, N_Drou_IDX
Character (Len=4) :: Drou_IDX(2), Dr_IDX
Logical :: status
!------------------------------------------------------------------------------
Drou_IDX(1)="PDSI"; Drou_IDX(2)="SPI"; N_Drou_IDX=2
!-----------------------------------------------------------------------------
!Dir="/home/tanpv/SEACLID_New/RES/POST"
!Dr_IDX="PDSI"
!Fname="/home/tanpv/SEACLID_New/Post_Pro_Stn_No_Rain_Center/PDSI.txt"

! Lay thong tin cac vung khi hau: B1=2, B2=3,..., N3=8
Call rd_ncdf_2d ("/home/tanpv/SEACLID/MASK/VN.grid_0.25.nc",loninfo,latinfo,nxlon,nxlat,vietnam,"all")
if (nxlon/=Nlon.OR.nxlat/=Nlat) then; print*, "Error in domain size !!!"; Stop; endif
Reg(1:Nlon,1:Nlat) = vietnam(1:Nlon,1:Nlat); Lon(1:Nlon) = loninfo(1:Nlon); Lat(1:Nlat) = latinfo(1:Nlat) 
Deallocate(loninfo, latinfo, vietnam)
Where (Reg>1) Reg=Reg-1

Open (1, file = "File_Drou_IDX.txt"); Read (1, Filenames); Close(1)    ! Doc ten file
Do iidx=1, 1! N_Drou_IDX
   Dr_IDX=trim(Drou_IDX(iidx))
   Do ifile=1, 1!NFile
      Fname=trim(Fname_inp(ifile))//trim(Dr_IDX)//".nc"
      if (iidx==1) then           ! PDSI
         Call Read_PDSI_nc (FName); PDSI(1:Nlon,1:Nlat,1:Ntime) = Dr_IDX_Tmp(1:Nlon,1:Nlat,1:Ntime); Deallocate(Dr_IDX_Tmp)
         Thre_1=Thre_1_PDSI; Thre_2=Thre_2_PDSI; Thre_3=Thre_3_PDSI; BIN(1:Nbin+1)=BIN_PDSI(1:Nbin+1)
      else                        ! SPI
         Call Read_SPI_nc (FName);  PDSI(1:Nlon,1:Nlat,1:Ntime) = Dr_IDX_Tmp(1:Nlon,1:Nlat,1:Ntime); Deallocate(Dr_IDX_Tmp)
         Thre_1=Thre_1_SPI; Thre_2=Thre_2_SPI; Thre_3=Thre_3_SPI; BIN(1:Nbin+1)=BIN_SPI(1:Nbin+1)
      endif
      Where (PDSI/=Miss.AND.PDSI<-20) PDSI=Miss    ! Avoid singular values
      Models=trim(Fname_Out(ifile))

      Call Cal_Freq_All_Timeseries(Dr_IDX, Models)      ! Tinh phan bo tan suat cho toan chuoi cua tung vung:    Freq(Nreg, Nbin)
      Call Cal_Freq_Monthly(Dr_IDX, Models)             ! Tinh phan bo tan suat cho tung thang cua tung vung:    Freq_Mon(Nreg, Nbin, Nmon)
      Call Drought_Index_Areal_Mean     ! Tinh trung binh chi so han cho tung vung:              PDSI_Reg(Ntime,Nreg)
      Call Define_Drought_Events        ! Xac dinh su kien han (0/1) tung thang cho tung tram va cho tung vung: D1, D2,...,OD1, OD2,..., D1_Reg,...,OD1_Reg,...
      Call Cal_Spells                   ! Tinh tan so theo do dai cac dot va tan so tung dot
!      Call Geographic_Extent            ! Tinh phan tram so tram bi han trong tung dot cua tung vung
      Call Geographic_Extent_1   ! Tinh phan tram so tram bi han trong tung THANG cua tung vung
      Call Cal_DS_DI                    ! Tinh Do khac nghiet (Severity) va Cuong do (Intensity) cua tung dot han
      Call Write_Spell_GeoExtent_To_NC_1(Dr_IDX, Models)  ! Write Spells & Geographic_Extent to NC file (xep theo thoi gian xuat hien)
      Call Write_Freq_Duration_To_NC    (Dr_IDX, Models)  ! Tan so cua cac dot theo vung va tung grid
      Call Write_Spell_To_NC(Dr_IDX, Models) ! Do dai cac dot cua tung vung va tung tram xep theo trinh tu
      Call Write_Severity_Intensity_To_NC(Dr_IDX, Models) ! Drought Severity & Intensity
   Enddo
Enddo
!Call Print_Test(Spel_Reg1,NSpel_Reg1)

CONTAINS
!-----------------------------------------------------------------------------
Subroutine Write_Spell_Area_To_TXT (X, Y, Fname)
Integer :: it, ireg, iyr, imon
Real, Dimension (Nreg,Ntime) :: X, Y
Character (Len=*) :: Fname
Open (3, file=trim(Fname),status="Unknown")
write (3, "(A8,7A16)") "YYYYMM", "R1","R2","R3","R4","R5","R6","R7"
write (3, "(15A8)") " ", "Spell","Area","Spell","Area","Spell","Area","Spell","Area","Spell","Area","Spell","Area","Spell","Area"
it=0
Do iyr=Yr1, Yr2
   Do imon=1, Nmon
      it=it+1
      write(3,"(1x,I4.4,I2.2,1x,20F8.3)") iyr, imon, (X(ireg,it),Y(ireg,it), ireg=1, Nreg) 
   Enddo
Enddo
Close(3)
End Subroutine Write_Spell_Area_To_TXT
!-----------------------------------------------------------------------------
Subroutine Write_Spell_GeoExtent_To_NC_1(Dr_IDX, Models)
Integer :: ireg, it, ispl, NSpel, NSpel_Max
Character (Len=*) :: Dr_IDX, Models
Real :: Xreg(Nreg), YSpel(Ntime), XArea(Nreg,Ntime), XSpell(Nreg,Ntime)  ! Thu tu vung, Thu tu cac dot, Geo-Extent, Duration of each spell
Real, Allocatable :: Xlon(:),Xlat(:), NumOfSpell(:,:), MaxSpell(:,:)   ! So dot han va do dai dai nhat cua cac dot han cua tung o luoi trong toan giai doan
!---------------------------------------
Allocate(Xlon(Nlon),Xlat(Nlat),NumOfSpell(Nlon,Nlat),MaxSpell(Nlon,Nlat))
Fname=trim(Models)//trim(Dr_IDX)//"_Num_Spell_and_Max_Leng_Spell_D1-3.nc"

Xlon(1:Nlon)=Lon(1:Nlon); Xlat(1:Nlat)=Lat(1:Nlat)
NumOfSpell(1:Nlon,1:Nlat)=N_Spel1(1:Nlon,1:Nlat); MaxSpell(1:Nlon,1:Nlat)=N_Dur1(1:Nlon,1:Nlat)
Call wr_ncdf_2d     (trim(Fname), Xlon, Xlat, Nlon, Nlat, NumOfSpell, "NSpl_D1", "Number of Drought Spells of D1", "Spell", "Lon", "Lat", Miss)
Call wr_ncdf_2d_apd (trim(Fname), MaxSpell, "MxLenSpl_D1", Nlon, Nlat, "Maximum length of Spells of D1", "Mon", "Lon", "Lat", Miss)

NumOfSpell(1:Nlon,1:Nlat)=N_Spel2(1:Nlon,1:Nlat); MaxSpell(1:Nlon,1:Nlat)=N_Dur2(1:Nlon,1:Nlat)
Call wr_ncdf_2d_apd (trim(Fname), NumOfSpell,  "NSpl_D2", Nlon, Nlat, "Number of Drought Spells of D2", "Mon", "Lon", "Lat", Miss)
Call wr_ncdf_2d_apd (trim(Fname), MaxSpell,"MxLenSpl_D2", Nlon, Nlat, "Maximum length of Spells of D2", "Mon", "Lon", "Lat", Miss)

NumOfSpell(1:Nlon,1:Nlat)=N_Spel3(1:Nlon,1:Nlat); MaxSpell(1:Nlon,1:Nlat)=N_Dur3(1:Nlon,1:Nlat)
Call wr_ncdf_2d_apd (trim(Fname), NumOfSpell,  "NSpl_D3", Nlon, Nlat, "Number of Drought Spells of D3", "Mon", "Lon", "Lat", Miss)
Call wr_ncdf_2d_apd (trim(Fname), MaxSpell,"MxLenSpl_D3", Nlon, Nlat, "Maximum length of Spells of D3", "Mon", "Lon", "Lat", Miss)

!---------------------------------------
Do ireg=1, Nreg; Xreg(ireg)=ireg; Enddo; Do it=1, Ntime; YSpel(it)=it; Enddo   ! Tao mang X, Y gia
! Drough D1
!XSpell=Miss; XArea=Miss
XSpell=0; XArea=0
Do ireg=1, Nreg
   NSpel=NSpel_Reg1(ireg)
   Do ispl=1, NSpel
      it=int((IS_Reg1(ireg,ispl)+IE_Reg1(ireg,ispl))/2.)
      XSpell(ireg,it) = Spel_Reg1(ireg,ispl); ! XArea(ireg,it) = Area_Reg1(ireg,ispl)
   Enddo
   XArea(ireg,1:Ntime) = Area_Reg1(ireg,1:Ntime)
Enddo
Fname=trim(Models)//trim(Dr_IDX)//"_Stat_Area_Spell_D1.nc"
Call wr_ncdf_2d (trim(Fname), Xreg, YSpel, Nreg, Ntime, XSpell, "Spells", "Drought Spells", "Mon", "Subregion", "Spell", Miss)
Call wr_ncdf_2d_apd (trim(Fname), XArea, "Area", Nreg, Ntime, "Geographic Extent", "%", "Subregion", "Spell", Miss)

!Fname=trim(Models)//trim(Dr_IDX)//"_Stat_Area_Spell_D1.txt"
!Call Write_Spell_Area_To_TXT (XSpell, XArea, Fname)

! Drough D2
XSpell=0; XArea=0 !XSpell=Miss; XArea=Miss
Do ireg=1, Nreg
   NSpel=NSpel_Reg2(ireg)
   Do ispl=1, NSpel
      it=int((IS_Reg2(ireg,ispl)+IE_Reg2(ireg,ispl))/2.)
      XSpell(ireg,it) = Spel_Reg2(ireg,ispl); ! XArea(ireg,it) = Area_Reg2(ireg,ispl)
   Enddo
   XArea(ireg,1:Ntime) = Area_Reg2(ireg,1:Ntime)
Enddo

Fname=trim(Models)//trim(Dr_IDX)//"_Stat_Area_Spell_D2.nc"
Call wr_ncdf_2d (trim(Fname), Xreg, YSpel, Nreg, Ntime, XSpell, "Spells", "Drought Spells", "Mon", "Subregion", "Spell", Miss)
Call wr_ncdf_2d_apd (trim(Fname), XArea, "Area", Nreg, Ntime, "Geographic Extent", "%", "Subregion", "Spell", Miss)

!Fname=trim(Models)//trim(Dr_IDX)//"_Stat_Area_Spell_D2.txt"
!Call Write_Spell_Area_To_TXT (XSpell, XArea, Fname)

! Drough D3
XSpell=0; XArea=0 !XSpell=Miss; XArea=Miss
Do ireg=1, Nreg
   NSpel=NSpel_Reg3(ireg)
   Do ispl=1, NSpel
      it=int((IS_Reg3(ireg,ispl)+IE_Reg3(ireg,ispl))/2.)
      XSpell(ireg,it) = Spel_Reg3(ireg,ispl);  !XArea(ireg,it) = Area_Reg3(ireg,ispl)
   Enddo
   XArea(ireg,1:Ntime) = Area_Reg3(ireg,1:Ntime)
Enddo

Fname=trim(Models)//trim(Dr_IDX)//"_Stat_Area_Spell_D3.nc"
Call wr_ncdf_2d (trim(Fname), Xreg, YSpel, Nreg, Ntime, XSpell, "Spells", "Drought Spells", "Mon", "Subregion", "Spell", Miss)
Call wr_ncdf_2d_apd (trim(Fname), XArea, "Area", Nreg, Ntime, "Geographic Extent", "%", "Subregion", "Spell", Miss)

!Fname=trim(Models)//trim(Dr_IDX)//"_Stat_Area_Spell_D3.txt"
!Call Write_Spell_Area_To_TXT (XSpell, XArea, Fname)

End Subroutine Write_Spell_GeoExtent_To_NC_1
!-----------------------------------------------------------------------------
Subroutine Geographic_Extent
Integer :: it, ilon, ilat, ireg, ispl, NSpel, ista, iend, idrou, itotal
! Area_Reg1(Nreg,Ntime)
 Area_Reg1=Miss; Area_Reg2=Miss; Area_Reg3=Miss
! Drough D1
Do ireg=1, Nreg
   NSpel=NSpel_Reg1(ireg)   ! Number of Spells
   Do ispl=1, NSpel
      idrou=0; itotal=0     ! Tong so tram bi han va Tong so tram duoc xet trong tung spell
      ista=IS_Reg1(ireg,ispl); iend=IE_Reg1(ireg,ispl)   ! Start and End time of this spell
      Do it=ista, iend      ! Search over entire spell
         Do ilon=1, Nlon
            Do ilat=1, Nlat
               if (Reg(ilon,ilat)==ireg) then
                  itotal=itotal+1        ! Tong so tram duoc xet trong tung spell
                  if (D1(ilon,ilat,it)==1) idrou=idrou+1   ! Tong so tram bi anh huong
               endif
            Enddo
         Enddo 
      Enddo
      Area_Reg1(ireg,ispl)=real(idrou)/real(itotal)*100.   ! Phan tram so tram bi anh huong
!      print*,ireg, Spel_Reg1(ireg,ispl), itotal, idrou, Area_Reg1(ireg,ispl); pause
   Enddo
Enddo
! Drough D2
Do ireg=1, Nreg
   NSpel=NSpel_Reg2(ireg)   ! Number of Spells
   Do ispl=1, NSpel
      idrou=0; itotal=0     ! Tong so tram bi han va Tong so tram duoc xet trong tung spell
      ista=IS_Reg2(ireg,ispl); iend=IE_Reg2(ireg,ispl)   ! Start and End time of this spell
      Do it=ista, iend      ! Search over entire spell
         Do ilon=1, Nlon
            Do ilat=1, Nlat
               if (Reg(ilon,ilat)==ireg) then
                  itotal=itotal+1        ! Tong so tram duoc xet trong tung spell
                  if (D2(ilon,ilat,it)==1) idrou=idrou+1   ! Tong so tram bi anh huong
               endif
            Enddo
         Enddo
      Enddo
      Area_Reg2(ireg,ispl)=real(idrou)/real(itotal)*100.   ! Phan tram so tram bi anh huong
!      print*,ireg, Spel_Reg1(ireg,ispl), itotal, idrou, Area_Reg1(ireg,ispl); pause
   Enddo
Enddo
! Drough D3
Do ireg=1, Nreg
   NSpel=NSpel_Reg3(ireg)   ! Number of Spells
   Do ispl=1, NSpel
      idrou=0; itotal=0     ! Tong so tram bi han va Tong so tram duoc xet trong tung spell
      ista=IS_Reg3(ireg,ispl); iend=IE_Reg3(ireg,ispl)   ! Start and End time of this spell
      Do it=ista, iend      ! Search over entire spell
         Do ilon=1, Nlon
            Do ilat=1, Nlat
               if (Reg(ilon,ilat)==ireg) then
                  itotal=itotal+1        ! Tong so tram duoc xet trong tung spell
                  if (D3(ilon,ilat,it)==1) idrou=idrou+1   ! Tong so tram bi anh huong
               endif
            Enddo
         Enddo
      Enddo
      Area_Reg3(ireg,ispl)=real(idrou)/real(itotal)*100.   ! Phan tram so tram bi anh huong
!      print*,ireg, Spel_Reg1(ireg,ispl), itotal, idrou, Area_Reg1(ireg,ispl); pause
   Enddo
Enddo
End Subroutine Geographic_Extent
!-----------------------------------------------------------------------------
Subroutine Print_Test(X, Y)
Integer :: i,j, N,M
Real :: X(Nreg,Ntime), Y(Nreg)
print*,"Testing"
Do i=1,Nreg
   print "(I4,100F4.0)", i, (X(i,j),j=1,Int(Y(i)))
Enddo
End Subroutine Print_Test
!-----------------------------------------------------------------------------
Subroutine Geographic_Extent_1
Integer :: it, ilon, ilat, ireg, idrou1,idrou2,idrou3, itotal
! Area_Reg1(Nreg,Ntime)
Area_Reg1=Miss; Area_Reg2=Miss; Area_Reg3=Miss
Do ireg=1, Nreg
   Do it=1, Ntime
      idrou1=0; idrou2=0; idrou3=0; itotal=0     ! Tong so tram bi han va Tong so tram duoc xet trong tung spell
      Do ilon=1, Nlon
         Do ilat=1, Nlat
            if (Int(Reg(ilon,ilat))==ireg) then
               itotal=itotal+1        ! Tong so tram duoc xet trong 1 thang  
               if (D1(ilon,ilat,it)==1) idrou1=idrou1+1   ! Tong so tram bi anh huong
               if (D2(ilon,ilat,it)==1) idrou2=idrou2+1   ! Tong so tram bi anh huong
               if (D3(ilon,ilat,it)==1) idrou3=idrou3+1   ! Tong so tram bi anh huong
            endif
         Enddo
      Enddo
      Area_Reg1(ireg,it)=real(idrou1)/real(itotal)*100.   ! Phan tram so tram bi anh huong
      Area_Reg2(ireg,it)=real(idrou2)/real(itotal)*100.   ! Phan tram so tram bi anh huong
      Area_Reg3(ireg,it)=real(idrou3)/real(itotal)*100.   ! Phan tram so tram bi anh huong
   Enddo
Enddo
End Subroutine Geographic_Extent_1
!-----------------------------------------------------------------------------
Subroutine Cal_Spells
Integer :: ilon, ilat, ireg, NDot, MaxDot
Real :: X(Ntime), M(Ntime), Spell(Ntime), I_Start(Ntime), I_End(Ntime), Imax, NSpell
!Integer :: Max_So_Dot        ! So dot lon nhat cua tat ca cac tram (=Max {Ndot} - Gia tri lon nhat cua N_Spell)
!Integer :: Max_Length_Dot    ! Do keo dai dai nhat trong tat ca cac dot cua cac tram (=Max {MaxDot} - Gia tri lon nhat cua Imax)

! M(MaxDot): Luu tan so cua tung dot, tu 1,2,..,MaxDot; MaxDot la do dai cuc dai cua cac dot
! Spell(NDot): Luu do dai tung dot, tu 1,2,..,NDot; NDot la so dot trong toan chuoi
! I_Start(NDot), I_End(NDot): Luu chi so dau va cuoi cua tung dot
Dur1=Miss; Dur2=Miss; Dur3=Miss;  Spel1=Miss; Spel2=Miss; Spel3=Miss; 
Max_So_Dot=0; Max_Length_Dot=0; Max_So_Dot_Reg=0; Max_Length_Dot_Reg=0
Do ilon=1, Nlon   ! Tinh tan so ung voi do dai tung dot va Tan so tung dot cua tung tram
   Do ilat=1, Nlat
      X(1:Ntime) = D1(ilon,ilat,1:Ntime); Call NSpell_and_Leng (X, Ntime, M, Imax, Spell, NSpell, I_Start, I_End); MaxDot=Int(Imax); NDot=Int(NSpell)
      Dur1(ilon,ilat,1:MaxDot)=M(1:MaxDot); Spel1(ilon,ilat,1:NDot)=Spell(1:NDot); IS1(ilon,ilat,1:NDot)=I_Start(1:NDot); IE1(ilon,ilat,1:NDot)=I_End(1:NDot)
      N_Dur1(ilon,ilat)=MaxDot; N_Spel1(ilon,ilat)=NDot
      if (Max_Length_Dot<MaxDot) Max_Length_Dot=MaxDot; if (Max_So_Dot<NDot) Max_So_Dot=NDot

      X(1:Ntime) = D2(ilon,ilat,1:Ntime); Call NSpell_and_Leng (X, Ntime, M, Imax, Spell, NSpell, I_Start, I_End); MaxDot=Int(Imax); NDot=Int(NSpell)
      Dur2(ilon,ilat,1:MaxDot)=M(1:MaxDot); Spel2(ilon,ilat,1:NDot)=Spell(1:NDot); IS2(ilon,ilat,1:NDot)=I_Start(1:NDot); IE2(ilon,ilat,1:NDot)=I_End(1:NDot)
      N_Dur2(ilon,ilat)=MaxDot; N_Spel2(ilon,ilat)=NDot
      if (Max_Length_Dot<MaxDot) Max_Length_Dot=MaxDot; if (Max_So_Dot<NDot) Max_So_Dot=NDot

      X(1:Ntime) = D3(ilon,ilat,1:Ntime); Call NSpell_and_Leng (X, Ntime, M, Imax, Spell, NSpell, I_Start, I_End); MaxDot=Int(Imax); NDot=Int(NSpell)
      Dur3(ilon,ilat,1:MaxDot)=M(1:MaxDot); Spel3(ilon,ilat,1:NDot)=Spell(1:NDot); IS3(ilon,ilat,1:NDot)=I_Start(1:NDot); IE3(ilon,ilat,1:NDot)=I_End(1:NDot)
      N_Dur3(ilon,ilat)=MaxDot; N_Spel3(ilon,ilat)=NDot
      if (Max_Length_Dot<MaxDot) Max_Length_Dot=MaxDot; if (Max_So_Dot<NDot) Max_So_Dot=NDot
   Enddo
Enddo
Dur_Reg1=Miss; Dur_Reg2=Miss; Dur_Reg3=Miss; Spel_Reg1=Miss; Spel_Reg2=Miss; Spel_Reg3=Miss
Do ireg=1, Nreg   ! Tinh tan so ung voi do dai tung dot va Tan so tung dot cua tung vung
   X(1:Ntime) = D1_Reg(ireg,1:Ntime); Call NSpell_and_Leng (X, Ntime, M, Imax, Spell, NSpell, I_Start, I_End); MaxDot=Int(Imax); NDot=Int(NSpell)
!print*,MaxDot,NDot, M(1:MaxDot), Spell(1:NDot); pause
   Dur_Reg1(ireg,1:MaxDot)=M(1:MaxDot); Spel_Reg1(ireg,1:NDot)=Spell(1:NDot); IS_Reg1(ireg,1:NDot)=I_Start(1:NDot); IE_Reg1(ireg,1:NDot)=I_End(1:NDot)
   NDur_Reg1(ireg)=MaxDot; NSpel_Reg1(ireg)=NDot
   if (Max_Length_Dot_Reg<MaxDot) Max_Length_Dot_Reg=MaxDot; if (Max_So_Dot_Reg<NDot) Max_So_Dot_Reg=NDot

   X(1:Ntime) = D2_Reg(ireg,1:Ntime); Call NSpell_and_Leng (X, Ntime, M, Imax, Spell, NSpell, I_Start, I_End); MaxDot=Int(Imax); NDot=Int(NSpell)
   Dur_Reg2(ireg,1:MaxDot)=M(1:MaxDot); Spel_Reg2(ireg,1:NDot)=Spell(1:NDot); IS_Reg2(ireg,1:NDot)=I_Start(1:NDot); IE_Reg2(ireg,1:NDot)=I_End(1:NDot)
   NDur_Reg2(ireg)=MaxDot; NSpel_Reg2(ireg)=NDot
   if (Max_Length_Dot_Reg<MaxDot) Max_Length_Dot_Reg=MaxDot; if (Max_So_Dot_Reg<NDot) Max_So_Dot_Reg=NDot

   X(1:Ntime) = D3_Reg(ireg,1:Ntime); Call NSpell_and_Leng (X, Ntime, M, Imax, Spell, NSpell, I_Start, I_End); MaxDot=Int(Imax); NDot=Int(NSpell)
   Dur_Reg3(ireg,1:MaxDot)=M(1:MaxDot); Spel_Reg3(ireg,1:NDot)=Spell(1:NDot); IS_Reg3(ireg,1:NDot)=I_Start(1:NDot); IE_Reg3(ireg,1:NDot)=I_End(1:NDot)
   NDur_Reg3(ireg)=MaxDot; NSpel_Reg3(ireg)=NDot
   if (Max_Length_Dot_Reg<MaxDot) Max_Length_Dot_Reg=MaxDot; if (Max_So_Dot_Reg<NDot) Max_So_Dot_Reg=NDot
Enddo
End Subroutine Cal_Spells
!-----------------------------------------------------------------------------
Subroutine Cal_DS_DI  ! Calculation of Drought Severity and Drought Intensity for each drought spell
Integer :: ilon, ilat, ireg, ispl, it, it1, it2, NDS1, NDS2, NDS3, NDI1, NDI2, NDI3
Real :: SumIDX, NSpl, Mean, Std
Real :: XDS1(Nlon*Nlat),XDS2(Nlon*Nlat),XDS3(Nlon*Nlat),XDI1(Nlon*Nlat),XDI2(Nlon*Nlat),XDI3(Nlon*Nlat)
! For grid cells
DS1=Miss; DS2=Miss; DS3=Miss; DI1=Miss; DI2=Miss; DI3=Miss
Do ilon=1, Nlon   ! Tinh tan so ung voi do dai tung dot va Tan so tung dot cua tung tram
   Do ilat=1, Nlat
      Do ispl=1, Int(N_Spel1(ilon,ilat))
         it1=IS1(ilon,ilat,ispl); it2=IE1(ilon,ilat,ispl); Nspl=Spel1(ilon,ilat,ispl); SumIDX=0.
         Do it=it1, it2
            if (PDSI(ilon,ilat,it)/=Miss) SumIDX=SumIDX+PDSI(ilon,ilat,it)     ! Tong tich luy chi so drought trong 1 spell
         Enddo
         DS1(ilon,ilat,ispl)=SumIDX; DI1(ilon,ilat,ispl)=SumIDX/Nspl    ! Drought Severity & Drought Intensity of each Spell
      Enddo
      Do ispl=1, Int(N_Spel2(ilon,ilat))
         it1=IS2(ilon,ilat,ispl); it2=IE2(ilon,ilat,ispl); Nspl=Spel2(ilon,ilat,ispl); SumIDX=0.
         Do it=it1, it2
            if (PDSI(ilon,ilat,it)/=Miss) SumIDX=SumIDX+PDSI(ilon,ilat,it)     ! Tong tich luy chi so drought trong 1 spell
         Enddo
         DS2(ilon,ilat,ispl)=SumIDX; DI2(ilon,ilat,ispl)=SumIDX/Nspl    ! Drought Severity & Drought Intensity of each Spell
      Enddo
      Do ispl=1, Int(N_Spel3(ilon,ilat))
         it1=IS3(ilon,ilat,ispl); it2=IE3(ilon,ilat,ispl); Nspl=Spel3(ilon,ilat,ispl); SumIDX=0.
         Do it=it1, it2
            if (PDSI(ilon,ilat,it)/=Miss) SumIDX=SumIDX+PDSI(ilon,ilat,it)     ! Tong tich luy chi so drought trong 1 spell
         Enddo
         DS3(ilon,ilat,ispl)=SumIDX; DI3(ilon,ilat,ispl)=SumIDX/Nspl    ! Drought Severity & Drought Intensity of each Spell
      Enddo
   Enddo
Enddo
! For subregions
! Ave_Std(X, N, Mean, Std)
DS1_Reg=Miss; DS2_Reg=Miss; DS2_Reg=Miss; DI1_Reg=Miss; DI2_Reg=Miss; DI3_Reg=Miss
DS1_Reg_Std=Miss; DS2_Reg_Std=Miss; DS3_Reg_Std=Miss; DI1_Reg_Std=Miss; DI2_Reg_Std=Miss; DI3_Reg_Std=Miss
Do ireg=1, Nreg
   Do ispl=1, Max_So_Dot
      NDS1=0; NDS2=0; NDS3=0; NDI1=0; NDI2=0; NDI3=0
      Do ilon=1, Nlon
         Do ilat=1, Nlat
            if (Int(Reg(ilon,ilat))==ireg) then
               if (DS1(ilon,ilat,ispl)/=Miss) then; NDS1=NDS1+1; XDS1(NDS1)=DS1(ilon,ilat,ispl); endif
               if (DS2(ilon,ilat,ispl)/=Miss) then; NDS2=NDS2+1; XDS2(NDS2)=DS2(ilon,ilat,ispl); endif
               if (DS3(ilon,ilat,ispl)/=Miss) then; NDS3=NDS3+1; XDS3(NDS3)=DS3(ilon,ilat,ispl); endif
               if (DI1(ilon,ilat,ispl)/=Miss) then; NDI1=NDI1+1; XDI1(NDI1)=DI1(ilon,ilat,ispl); endif
               if (DI2(ilon,ilat,ispl)/=Miss) then; NDI2=NDI2+1; XDI2(NDI2)=DI2(ilon,ilat,ispl); endif
               if (DI3(ilon,ilat,ispl)/=Miss) then; NDI3=NDI3+1; XDI3(NDI3)=DI3(ilon,ilat,ispl); endif
            endif
         Enddo
      Enddo
      Call Ave_Std(XDS1, NDS1, Mean, Std); DS1_Reg(ireg,ispl)=Mean; DS1_Reg_Std(ireg,ispl)=Std
      Call Ave_Std(XDS2, NDS2, Mean, Std); DS2_Reg(ireg,ispl)=Mean; DS2_Reg_Std(ireg,ispl)=Std
      Call Ave_Std(XDS3, NDS3, Mean, Std); DS3_Reg(ireg,ispl)=Mean; DS3_Reg_Std(ireg,ispl)=Std
      Call Ave_Std(XDI1, NDI1, Mean, Std); DI1_Reg(ireg,ispl)=Mean; DI1_Reg_Std(ireg,ispl)=Std
      Call Ave_Std(XDI2, NDI2, Mean, Std); DI2_Reg(ireg,ispl)=Mean; DI2_Reg_Std(ireg,ispl)=Std
      Call Ave_Std(XDI3, NDI3, Mean, Std); DI3_Reg(ireg,ispl)=Mean; DI3_Reg_Std(ireg,ispl)=Std
   Enddo
Enddo

!DS1_Reg=Miss; DS2_Reg=Miss; DS2_Reg=Miss; DI1_Reg=Miss; DI2_Reg=Miss; DI3_Reg=Miss
!Do ireg=1, Nreg 
!   Do ispl=1, Int(NSpel_Reg1(ireg))
!      it1=IS_Reg1(ireg,ispl); it2=IE_Reg1(ireg,ispl); Nspl=Spel_Reg1(ireg,ispl); SumIDX=0.
!      Do it=it1, it2
!         if (PDSI_Reg(it,ireg)/=Miss) SumIDX=SumIDX+PDSI_Reg(it,ireg)      ! Tong tich luy chi so drought trong 1 spell
!      Enddo
!      DS1_Reg(ireg,ispl)=SumIDX; DI1_Reg(ireg,ispl)=SumIDX/Nspl  ! Drought Severity & Drought Intensity of each Spell
!   Enddo
!   Do ispl=1, Int(NSpel_Reg2(ireg))
!      it1=IS_Reg2(ireg,ispl); it2=IE_Reg2(ireg,ispl); Nspl=Spel_Reg2(ireg,ispl); SumIDX=0.
!      Do it=it1, it2
!         if (PDSI_Reg(it,ireg)/=Miss) SumIDX=SumIDX+PDSI_Reg(it,ireg)      ! Tong tich luy chi so drought trong 1 spell
!      Enddo
!      DS2_Reg(ireg,ispl)=SumIDX; DI2_Reg(ireg,ispl)=SumIDX/Nspl  ! Drought Severity & Drought Intensity of each Spell
!   Enddo
!   Do ispl=1, Int(NSpel_Reg3(ireg))
!      it1=IS_Reg3(ireg,ispl); it2=IE_Reg3(ireg,ispl); Nspl=Spel_Reg3(ireg,ispl); SumIDX=0.
!      Do it=it1, it2
!         if (PDSI_Reg(it,ireg)/=Miss) SumIDX=SumIDX+PDSI_Reg(it,ireg)      ! Tong tich luy chi so drought trong 1 spell
!      Enddo
!      DS3_Reg(ireg,ispl)=SumIDX; DI3_Reg(ireg,ispl)=SumIDX/Nspl  ! Drought Severity & Drought Intensity of each Spell
!   Enddo
!Enddo
End Subroutine Cal_DS_DI
!-----------------------------------------------------------------------------
Subroutine Write_Severity_Intensity_To_NC(Dr_IDX, Models)
Character (Len=*) :: Dr_IDX, Models
Real, Allocatable :: XDot(:)      ! So thu tu cac dot (=1,2,... Max_So_Dot)
Real, Allocatable :: Seve(:,:), Inten(:,:) ! Severity & Intensity cua tung vung (1:Nreg, 1:Max_So_Dot_Reg)
Real, Allocatable :: SVstd(:,:), INstd(:,:) ! Std of Severity & Intensity cua tung vung (1:Nreg, 1:Max_So_Dot_Reg)
Real, Allocatable :: DS_grd(:,:,:),DI_grd(:,:,:) ! Severity & Intensity cua tung grid (Nlon,Nlat,Max_So_Dot)
Real :: Xreg(Nreg)
Integer :: ilon,ilat, ireg, idot, LengDot
! For Sub-Regions:
!Allocate(XDot(1:Max_So_Dot_Reg),Seve(1:Nreg,1:Max_So_Dot_Reg),Inten(1:Nreg,1:Max_So_Dot_Reg))
Allocate(XDot(1:Max_So_Dot),Seve(1:Nreg,1:Max_So_Dot),Inten(1:Nreg,1:Max_So_Dot),SVstd(1:Nreg,1:Max_So_Dot),INstd(1:Nreg,1:Max_So_Dot))
Fname=trim(Models)//trim(Dr_IDX)//"_Severity_Intensity_Reg_D1-3.nc"

Do ireg=1, Nreg; Xreg(ireg)=ireg; Enddo                 ! Mang gia X Cac vung
Do idot=1, Max_So_Dot; XDot(idot)=idot; Enddo   ! Mang gia Y Do dai cac dot
! Drought D1:
Seve(:,:)=Miss; Inten(:,:)=Miss; SVstd(:,:)=Miss; INstd(:,:)=Miss; LengDot=Max_So_Dot
Do ireg=1, Nreg
!   LengDot=NSpel_Reg1(ireg); 
   Seve(ireg,1:LengDot)=DS1_Reg(ireg,1:LengDot); Inten(ireg,1:LengDot)=DI1_Reg(ireg,1:LengDot)
   SVstd(ireg,1:LengDot)=DS1_Reg_Std(ireg,1:LengDot); INstd(ireg,1:LengDot)=DI1_Reg_Std(ireg,1:LengDot)
Enddo
Call wr_ncdf_2d     (trim(Fname), Xreg, XDot, Nreg, Max_So_Dot, Seve, "Seve_D1", "Drought Severity of D1","-","Subregion","Spell", Miss)
Call wr_ncdf_2d_apd (trim(Fname), SVstd,"SVstd_D1",Nreg, Max_So_Dot,"Std of Dr. Severity of D1","-","Subregion","Spell", Miss)
Call wr_ncdf_2d_apd (trim(Fname), Inten,"Inte_D1", Nreg, Max_So_Dot,"Drought Intensity of D1", "-","Subregion","Spell", Miss)
Call wr_ncdf_2d_apd (trim(Fname), INstd,"INstd_D1",Nreg, Max_So_Dot,"Std of Dr. Intensity of D1","-","Subregion","Spell", Miss)
! Drought D2:
Seve(:,:)=Miss; Inten(:,:)=Miss; SVstd(:,:)=Miss; INstd(:,:)=Miss;
Do ireg=1, Nreg
!   LengDot=NSpel_Reg2(ireg); Seve(ireg,1:LengDot)=DS2_Reg(ireg,1:LengDot); Inten(ireg,1:LengDot)=DI2_Reg(ireg,1:LengDot)
   Seve(ireg,1:LengDot)=DS2_Reg(ireg,1:LengDot); Inten(ireg,1:LengDot)=DI2_Reg(ireg,1:LengDot)
   SVstd(ireg,1:LengDot)=DS2_Reg_Std(ireg,1:LengDot); INstd(ireg,1:LengDot)=DI2_Reg_Std(ireg,1:LengDot)
Enddo
Call wr_ncdf_2d_apd (trim(Fname), Seve, "Seve_D2",Nreg, Max_So_Dot, "Drought Severity of D2","-","Subregion","Spell", Miss)
Call wr_ncdf_2d_apd (trim(Fname), SVstd,"SVstd_D2",Nreg, Max_So_Dot,"Std of Dr. Severity of D2","-","Subregion","Spell", Miss)
Call wr_ncdf_2d_apd (trim(Fname), Inten,"Inte_D2",Nreg, Max_So_Dot,"Drought Intensity of D2","-","Subregion","Spell", Miss)
Call wr_ncdf_2d_apd (trim(Fname), INstd,"INstd_D2",Nreg, Max_So_Dot,"Std of Dr. Intensity of D2","-","Subregion","Spell", Miss)
! Drought D3:
Seve(:,:)=Miss; Inten(:,:)=Miss; SVstd(:,:)=Miss; INstd(:,:)=Miss;
Do ireg=1, Nreg
!   LengDot=NSpel_Reg3(ireg); Seve(ireg,1:LengDot)=DS3_Reg(ireg,1:LengDot); Inten(ireg,1:LengDot)=DI3_Reg(ireg,1:LengDot)
   Seve(ireg,1:LengDot)=DS3_Reg(ireg,1:LengDot); Inten(ireg,1:LengDot)=DI3_Reg(ireg,1:LengDot)
   SVstd(ireg,1:LengDot)=DS3_Reg_Std(ireg,1:LengDot); INstd(ireg,1:LengDot)=DI3_Reg_Std(ireg,1:LengDot)
Enddo
Call wr_ncdf_2d_apd (trim(Fname), Seve, "Seve_D3",Nreg, Max_So_Dot, "Drought Severity of D3","-","Subregion","Spell", Miss)
Call wr_ncdf_2d_apd (trim(Fname), SVstd,"SVstd_D3",Nreg, Max_So_Dot,"Std of Dr. Severity of D3","-","Subregion","Spell", Miss)
Call wr_ncdf_2d_apd (trim(Fname), Inten,"Inte_D3",Nreg, Max_So_Dot,"Drought Intensity of D3","-","Subregion","Spell", Miss)
Call wr_ncdf_2d_apd (trim(Fname), INstd,"INstd_D3",Nreg, Max_So_Dot,"Std of Dr. Intensity of D3","-","Subregion","Spell", Miss)

! For each Station:
Fname=trim(Models)//trim(Dr_IDX)//"_Severity_Intensity_Grd_D1-3.nc"
DeAllocate(XDot,Seve,Inten); Allocate(DS_grd(Nlon,Nlat,Max_So_Dot),DI_grd(Nlon,Nlat,Max_So_Dot),XDot(Max_So_Dot))  ! Mang gia XDot
Do idot=1, Max_So_Dot; XDot(idot)=idot; Enddo

! Drought D1: 
DS_grd(:,:,:)=Miss; DI_grd(:,:,:)=Miss
Do ilon=1, Nlon
   Do ilat=1, Nlat
      LengDot=N_Spel1(ilon,ilat); DS_grd(ilon,ilat,1:LengDot)=DS1(ilon,ilat,1:LengDot); DI_grd(ilon,ilat,1:LengDot)=DI1(ilon,ilat,1:LengDot)
   Enddo
Enddo
Call wr_ncdf_3d(trim(Fname),Lon,Lat,XDot, Nlon,Nlat,Max_So_Dot, DS_grd,"Seve_D1","months since 1980-01-15 00:00:00","Severity of D1","-","Lon","Lat","Spell",Miss)
Call wr_ncdf_3d_apd(trim(Fname),DI_grd,"Inte_D1", Nlon,Nlat,Max_So_Dot, "Intensity of D1","-","Lon","Lat","Spell",Miss)
! Drought D2: 
DS_grd(:,:,:)=Miss; DI_grd(:,:,:)=Miss
Do ilon=1, Nlon
   Do ilat=1, Nlat
      LengDot=N_Spel2(ilon,ilat); DS_grd(ilon,ilat,1:LengDot)=DS2(ilon,ilat,1:LengDot); DI_grd(ilon,ilat,1:LengDot)=DI2(ilon,ilat,1:LengDot)
   Enddo
Enddo
Call wr_ncdf_3d_apd(trim(Fname),DS_grd,"Seve_D2", Nlon,Nlat,Max_So_Dot, "Severity of D2","-","Lon","Lat","Spell",Miss)
Call wr_ncdf_3d_apd(trim(Fname),DI_grd,"Inte_D2", Nlon,Nlat,Max_So_Dot,"Intensity of D2","-","Lon","Lat","Spell",Miss)
! Drought D2: 
DS_grd(:,:,:)=Miss; DI_grd(:,:,:)=Miss
Do ilon=1, Nlon
   Do ilat=1, Nlat
      LengDot=N_Spel3(ilon,ilat); DS_grd(ilon,ilat,1:LengDot)=DS3(ilon,ilat,1:LengDot); DI_grd(ilon,ilat,1:LengDot)=DI3(ilon,ilat,1:LengDot)
   Enddo
Enddo
Call wr_ncdf_3d_apd(trim(Fname),DS_grd,"Seve_D3", Nlon,Nlat,Max_So_Dot, "Severity of D3","-","Lon","Lat","Spell",Miss)
Call wr_ncdf_3d_apd(trim(Fname),DI_grd,"Inte_D3", Nlon,Nlat,Max_So_Dot,"Intensity of D3","-","Lon","Lat","Spell",Miss)

End Subroutine Write_Severity_Intensity_To_NC
!-----------------------------------------------------------------------------
Subroutine Write_Spell_To_NC(Dr_IDX, Models)
! Write Spells (Duration of each Spell)
Character (Len=*) :: Dr_IDX, Models
Real, Allocatable :: XDot(:)      ! So thu tu cac dot (=1,2,... Max_So_Dot)
Real, Allocatable :: Fre_Dur(:,:) ! Do dai cac dot cua tung vung (1:Nreg, 1:Max_So_Dot_Reg)
Real, Allocatable :: FD_grd(:,:,:) ! Do dai cac dot cua tung grid (Nlon,Nlat,Max_So_Dot)
Real :: Xreg(Nreg)
Integer :: ilon,ilat, ireg, idot, LengDot

! For Sub-Regions:
Allocate(XDot(1:Max_So_Dot_Reg),Fre_Dur(1:Nreg,1:Max_So_Dot_Reg))
Fname=trim(Models)//trim(Dr_IDX)//"_Spell_Duration_Reg_D1-3.nc"
! Drought D1:
Do ireg=1, Nreg; Xreg(ireg)=ireg; Enddo                 ! Mang gia X Cac vung
Do idot=1, Max_So_Dot_Reg; XDot(idot)=idot; Enddo   ! Mang gia Y Do dai cac dot
Fre_Dur(:,:)=Miss;
Do ireg=1, Nreg
   LengDot=NSpel_Reg1(ireg)
   Fre_Dur(ireg,1:LengDot)=Spel_Reg1(ireg,1:LengDot)
Enddo
Call wr_ncdf_2d (trim(Fname), Xreg, XDot, Nreg, Max_So_Dot_Reg, Fre_Dur, "Spell_D1", "Spell Duration of D1","Mon","Subregion","Spell", Miss)

! Drought D2:
Fre_Dur(:,:)=Miss;
Do ireg=1, Nreg
   LengDot=NSpel_Reg2(ireg)
   Fre_Dur(ireg,1:LengDot)=Spel_Reg2(ireg,1:LengDot)
Enddo
Call wr_ncdf_2d_apd (trim(Fname), Fre_Dur,"Spell_D2",Nreg, Max_So_Dot_Reg,"Spell Duration of D2","Mon","Subregion","Spell", Miss)

! Drought D3:
Fre_Dur(:,:)=Miss;
Do ireg=1, Nreg
   LengDot=NSpel_Reg3(ireg)
   Fre_Dur(ireg,1:LengDot)=Spel_Reg3(ireg,1:LengDot)
Enddo
Call wr_ncdf_2d_apd (trim(Fname), Fre_Dur,"Spell_D3",Nreg, Max_So_Dot_Reg,"Spell Duration of D3","Mon","Subregion","Spell", Miss)

! For each Station:

Fname=trim(Models)//trim(Dr_IDX)//"_Spell_Duration_Grd_D1-3.nc"
DeAllocate(XDot,Fre_Dur); Allocate(FD_grd(Nlon,Nlat,Max_So_Dot),XDot(Max_So_Dot))  ! Mang gia XDot
Do idot=1, Max_So_Dot; XDot(idot)=idot; Enddo

! Drought D1:
FD_grd(:,:,:)=Miss
Do ilon=1, Nlon
   Do ilat=1, Nlat
      LengDot=N_Spel1(ilon,ilat); FD_grd(ilon,ilat,1:LengDot)=Spel1(ilon,ilat,1:LengDot)
   Enddo
Enddo
!wr_ncdf_3d (filenm,lon,lat,t, nlon,nlat,nt, var,var_nm, t_unit, var_long_name, var_unit, lon_nm, lat_nm, t_nm, flag)
!wr_ncdf_3d_apd (filenm,var,var_nm,nlon,nlat,nt, var_long_name, var_unit, lon_nm, lat_nm, t_nm, flag)
Call wr_ncdf_3d(trim(Fname),Lon,Lat,XDot, Nlon,Nlat,Max_So_Dot, FD_grd,"Spell_D1","months since 1980-01-15 00:00:00","Spell Duration of D1","Mon","Lon","Lat","Spell",Miss)

! Drought D2:
FD_grd(:,:,:)=Miss
Do ilon=1, Nlon
   Do ilat=1, Nlat
      LengDot=N_Spel2(ilon,ilat); FD_grd(ilon,ilat,1:LengDot)=Spel2(ilon,ilat,1:LengDot)
   Enddo
Enddo
Call wr_ncdf_3d_apd(trim(Fname),FD_grd,"Spell_D2", Nlon,Nlat,Max_So_Dot, "Spell Duration of D2","Mon","Lon","Lat","Spell",Miss)

! Drought D3:
FD_grd(:,:,:)=Miss
Do ilon=1, Nlon
   Do ilat=1, Nlat
      LengDot=N_Spel3(ilon,ilat); FD_grd(ilon,ilat,1:LengDot)=Spel3(ilon,ilat,1:LengDot)
   Enddo
Enddo
Call wr_ncdf_3d_apd(trim(Fname),FD_grd,"Freq_D3", Nlon,Nlat,Max_So_Dot, "Spell Duration of D3","Mon","Lon","Lat","Spell",Miss)

End Subroutine Write_Spell_To_NC
!-----------------------------------------------------------------------------
Subroutine Write_Freq_Duration_To_NC(Dr_IDX, Models)
! Write Frequency of Duration (Tan so cua cac dot keo dai)
Character (Len=*) :: Dr_IDX, Models
Real, Allocatable :: XDot(:)      ! Do keo dai cac dot (=1,2,... Max_Length_Dot) (month)
Real, Allocatable :: Fre_Dur(:,:) ! Tan so cac dot cua tung vung (1:Nreg, 1:Max_Length_Dot)
Real, Allocatable :: FD_grd(:,:,:) ! Tan so cac dot cua tung grid (Nlon,Nlat,Max_Length_Dot)
Real :: Xreg(Nreg)
Integer :: ilon,ilat, ireg, idot, LengDot
! For Sub-Regions:
Allocate(XDot(1:Max_Length_Dot_Reg),Fre_Dur(1:Nreg,1:Max_Length_Dot_Reg)) 
Fname=trim(Models)//trim(Dr_IDX)//"_Freq_Duration_Reg_D1-3.nc"
! Drought D1:
Do ireg=1, Nreg; Xreg(ireg)=ireg; Enddo                 ! Mang gia X Cac vung
Do idot=1, Max_Length_Dot_Reg; XDot(idot)=idot; Enddo   ! Mang gia Y Do dai cac dot
Fre_Dur(:,:)=Miss; 
Do ireg=1, Nreg
   LengDot=NDur_Reg1(ireg)
   Fre_Dur(ireg,1:LengDot)=Dur_Reg1(ireg,1:LengDot)
Enddo
Call wr_ncdf_2d (trim(Fname), Xreg, XDot, Nreg, Max_Length_Dot_Reg, Fre_Dur, "Freq_D1", "Freq. of Duration of D1","So_Dot","Subregion","Month", Miss)

! Drought D2:
Fre_Dur(:,:)=Miss;
Do ireg=1, Nreg
   LengDot=NDur_Reg2(ireg)
   Fre_Dur(ireg,1:LengDot)=Dur_Reg2(ireg,1:LengDot)
Enddo
Call wr_ncdf_2d_apd (trim(Fname), Fre_Dur,"Freq_D2",Nreg, Max_Length_Dot_Reg,"Freq. of Duration of D2","So_Dot","Subregion","Month", Miss)

! Drought D3:
Fre_Dur(:,:)=Miss;
Do ireg=1, Nreg
   LengDot=NDur_Reg3(ireg)
   Fre_Dur(ireg,1:LengDot)=Dur_Reg3(ireg,1:LengDot)
Enddo
Call wr_ncdf_2d_apd (trim(Fname), Fre_Dur,"Freq_D3",Nreg, Max_Length_Dot_Reg,"Freq. of Duration of D3","So_Dot","Subregion","Month", Miss)

! For each Station:
Fname=trim(Models)//trim(Dr_IDX)//"_Freq_Duration_Grd_D1-3.nc"
DeAllocate(XDot,Fre_Dur); Allocate(FD_grd(Nlon,Nlat,Max_Length_Dot),XDot(Max_Length_Dot))  ! Mang gia XDot
Do idot=1, Max_Length_Dot; XDot(idot)=idot; Enddo

! Drought D1:
FD_grd(:,:,:)=Miss
Do ilon=1, Nlon
   Do ilat=1, Nlat
      LengDot=N_Dur1(ilon,ilat); FD_grd(ilon,ilat,1:LengDot)=Dur1(ilon,ilat,1:LengDot)
   Enddo
Enddo
!wr_ncdf_3d (filenm,lon,lat,t, nlon,nlat,nt, var,var_nm, t_unit, var_long_name, var_unit, lon_nm, lat_nm, t_nm, flag)
!wr_ncdf_3d_apd (filenm,var,var_nm,nlon,nlat,nt, var_long_name, var_unit, lon_nm, lat_nm, t_nm, flag)
Call wr_ncdf_3d(trim(Fname),Lon,Lat,XDot, Nlon,Nlat,Max_Length_Dot, FD_grd,"Freq_D1","months since 1980-01-15 00:00:00","Freq. of Duration of D1","Dot","Lon","Lat","Duration",Miss)

! Drought D2:
FD_grd(:,:,:)=Miss
Do ilon=1, Nlon
   Do ilat=1, Nlat
      LengDot=N_Dur2(ilon,ilat); FD_grd(ilon,ilat,1:LengDot)=Dur2(ilon,ilat,1:LengDot)
   Enddo
Enddo
Call wr_ncdf_3d_apd(trim(Fname),FD_grd,"Freq_D2", Nlon,Nlat,Max_Length_Dot, "Freq. of Duration of D2","Dot","Lon","Lat","Duration",Miss)

! Drought D3:
FD_grd(:,:,:)=Miss
Do ilon=1, Nlon
   Do ilat=1, Nlat
      LengDot=N_Dur3(ilon,ilat); FD_grd(ilon,ilat,1:LengDot)=Dur3(ilon,ilat,1:LengDot)
   Enddo
Enddo
Call wr_ncdf_3d_apd(trim(Fname),FD_grd,"Freq_D3", Nlon,Nlat,Max_Length_Dot, "Freq. of Duration of D3","Dot","Lon","Lat","Duration",Miss)

End Subroutine Write_Freq_Duration_To_NC
!-----------------------------------------------------------------------------
Subroutine Define_Drought_Events
! D1: Mild or more intense drought, D2: Moderate more intense drought, D3: Severe of more intense drought, D4: Extreme drought
Integer :: ilon, ilat, ireg, it, DR1, DR0
D1=Miss; D2=Miss; D3=Miss; OD1=Miss; OD2=Miss; OD3=Miss
Do ilon=1, Nlon
   Do ilat=1, Nlat
      Do it=1, Ntime
         if (PDSI(ilon,ilat,it)/=Miss.AND.PDSI(ilon,ilat,it)<=Thre_1) then; D1(ilon,ilat,it)=1; OD1(ilon,ilat,it)=0; else; D1(ilon,ilat,it)=0; OD1(ilon,ilat,it)=1; endif
         if (PDSI(ilon,ilat,it)/=Miss.AND.PDSI(ilon,ilat,it)<=Thre_2) then; D2(ilon,ilat,it)=1; OD2(ilon,ilat,it)=0; else; D2(ilon,ilat,it)=0; OD2(ilon,ilat,it)=1; endif
         if (PDSI(ilon,ilat,it)/=Miss.AND.PDSI(ilon,ilat,it)<=Thre_3) then; D3(ilon,ilat,it)=1; OD3(ilon,ilat,it)=0; else; D3(ilon,ilat,it)=0; OD3(ilon,ilat,it)=1; endif
      Enddo
   Enddo
Enddo
D1_reg=Miss; D2_reg=Miss; D3_reg=Miss; OD1_reg=Miss; OD2_reg=Miss; OD3_reg=Miss
Do ireg=1, Nreg
   Do it=1, Ntime
      Call Reg_Drought_Event(PDSI,ireg,it,Thre_1, Dr1, Dr0); D1_Reg(ireg,it)=Dr1; OD1_Reg(ireg,it)=Dr0
      Call Reg_Drought_Event(PDSI,ireg,it,Thre_2, Dr1, Dr0); D2_Reg(ireg,it)=Dr1; OD2_Reg(ireg,it)=Dr0
      Call Reg_Drought_Event(PDSI,ireg,it,Thre_3, Dr1, Dr0); D3_Reg(ireg,it)=Dr1; OD3_Reg(ireg,it)=Dr0
!      if (PDSI_Reg(it,ireg)/=Miss.AND.PDSI_Reg(it,ireg)<=Thre_1) then; D1_Reg(ireg,it)=1; OD1_Reg(ireg,it)=0; else; D1_Reg(ireg,it)=0; OD1_Reg(ireg,it)=1; endif
!      if (PDSI_Reg(it,ireg)/=Miss.AND.PDSI_Reg(it,ireg)<=Thre_2) then; D2_Reg(ireg,it)=1; OD2_Reg(ireg,it)=0; else; D2_Reg(ireg,it)=0; OD2_Reg(ireg,it)=1; endif
!      if (PDSI_Reg(it,ireg)/=Miss.AND.PDSI_Reg(it,ireg)<=Thre_3) then; D3_Reg(ireg,it)=1; OD3_Reg(ireg,it)=0; else; D3_Reg(ireg,it)=0; OD3_Reg(ireg,it)=1; endif
   Enddo
Enddo
End Subroutine Define_Drought_Events
!-----------------------------------------------------------------------------
Subroutine Drought_Index_Areal_Mean
Integer :: ireg, it
Do it=1, Ntime
   Do ireg=1, Nreg
      PDSI_Reg(it,ireg)=Area_Mean (PDSI, ireg, it)
   Enddo
Enddo
End Subroutine Drought_Index_Areal_Mean
!-----------------------------------------------------------------------------
Subroutine Cal_Freq_All_Timeseries(Dr_IDX, Models)
Integer :: it, ilon, ilat, ireg, ibin, ic, N
Real :: tmp
Character (Len=*) :: Dr_IDX, Models
Do ireg=1, Nreg
   Do ibin=1, Nbin
      ic=0; N=0
      Do ilon=1, Nlon
         Do ilat=1, Nlat
            if (Int(Reg(ilon,ilat))==ireg) then
               Do it=1, Ntime
                  if (PDSI(ilon,ilat,it) /= Miss) then
                     N=N+1
                     if (PDSI(ilon,ilat,it)>=BIN(ibin) .AND. PDSI(ilon,ilat,it)<BIN(ibin+1)) ic=ic+1
                  endif
               Enddo
            endif
         Enddo
      Enddo
      Freq(ireg,ibin) = Real(ic)/Real(N)*100; Xbin(ibin)=(BIN(ibin)+BIN(ibin+1))/2.
   !   print "(I2,2F6.2)",ireg, Xbin(ibin), Freq(ireg,ibin)
   Enddo
Enddo
! write to netCDF file:
!Fname="/home/tanpv/SEACLID_New/Post_Pro_Stn_No_Rain_Center/PDSI_Stat.nc"
Fname=trim(Models)//trim(Dr_IDX)//"_Reg_Fre_Distribution.nc"
Do ireg=1, Nreg; Xreg(ireg)=ireg; Enddo
Call wr_ncdf_2d (trim(Fname), Xreg, Xbin, Nreg, Nbin, Freq, "Freq", "Freq_Drought_Index", "-", "Sub-Region", "BIN", Miss)

End Subroutine Cal_Freq_All_Timeseries
!-----------------------------------------------------------------------------
Subroutine Cal_Freq_Monthly(Dr_IDX, Models)
Integer :: imon, imon1, it, ilon, ilat, ireg, ibin, ic, N
Real :: Xmon(Nmon)
Character (Len=*) :: Dr_IDX, Models
Character (Len=3), Parameter :: Mon_Name(Nmon)=(/"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"/)
Do ireg=1, Nreg
   Do ibin=1, Nbin
      Do imon=1, Nmon
         ic=0; N=0
         Do ilon=1, Nlon
            Do ilat=1, Nlat
               if (Int(Reg(ilon,ilat))==ireg) then
                  Do it=1, Ntime
                     if (PDSI(ilon,ilat,it) /= Miss) then
                        imon1=mod(it, Nmon); if (imon1==0) imon1=12 
                        if (imon1==imon) then
                           N=N+1
                           if (PDSI(ilon,ilat,it)>=BIN(ibin) .AND. PDSI(ilon,ilat,it)<BIN(ibin+1)) ic=ic+1
                        endif
                     endif
                  Enddo
               endif
            Enddo
         Enddo
         if (N>0) then; Freq_Mon(ireg,ibin,imon) = Real(ic)/Real(N)*100; else; Freq_Mon(ireg,ibin,imon) = Miss; endif
      !   print "(3I3,2F6.2)",ireg, ibin, imon, Freq_Mon(ireg,ibin,imon)
      Enddo
   Enddo
!   print*,Freq_Mon(ireg,1:Nbin,1:Nmon); pause
Enddo
Do ibin=1, Nbin; Xbin(ibin)=(BIN(ibin)+BIN(ibin+1))/2.; Enddo; Do imon=1, Nmon; Xmon(imon)=imon; Enddo
! write to netCDF file:

Fname=trim(Models)//trim(Dr_IDX)//"_Monthly_Reg_Fre_Distribution.nc"
Do ireg=1, Nreg; Xreg(ireg)=ireg; Enddo
!Call wr_ncdf_2d (trim(Fname), Xreg, Xbin, Nreg, Nbin, Freq, "Freq", "Freq_Drought_Index", "-", "Sub-Region", "BIN", Miss)
! wr_ncdf_3d (filenm,lon,lat,t,nlon,nlat,nt, var,var_nm,t_unit, var_long_name, var_unit, lon_nm, lat_nm, t_nm, flag)
Call wr_ncdf_3d (trim(FName), Xreg, Xbin, Xmon, Nreg, Nbin, Nmon, Freq_Mon, "Freq","months since 1980-01-15 00:00:00", &
                              "Mon Freq Distri", "%", "Subregion","BIN","Mon", Miss)
End Subroutine Cal_Freq_Monthly
!--------------------------------------------------------------------------------
Subroutine Reg_Drought_Event(X, ireg, it, Thre, Drought, No_Drought)
! Ham xac dinh thang thu it co Drought Event hay khong. Neu 50% so tram co Drought thi Drought=1, No_Drought=0
! Drought=1/0 and No_Drought=0/1
Integer :: it, ilon, ilat, ireg, idr, igrd, Drought, No_Drought
Real :: Thre, Area_Percent
Real, Dimension(Nlon,Nlat,Ntime) :: X
idr=0; igrd=0
Do ilon=1, Nlon
   Do ilat=1, Nlat
      if (INT(Reg(ilon,ilat))==ireg) then
         igrd=igrd+1
         if (X(ilon,ilat,it)/=Miss.AND.X(ilon,ilat,it)<Thre) idr=idr+1
      endif
   Enddo
Enddo
Area_Percent=Real(idr)/Real(igrd)
if (Area_Percent>0.5) then; Drought=1; No_Drought=0; else; Drought=0; No_Drought=1; endif
End Subroutine Reg_Drought_Event
!--------------------------------------------------------------------------------
Real Function Area_Mean (X, ireg, it)   ! Tinh trung binh vung chi so X nao do
Integer :: it, ilon, ilat, ireg, icount
Real, Dimension(Nlon,Nlat,Ntime) :: X
Real :: tmp
tmp = 0; icount=0
Do ilon=1, Nlon
   Do ilat=1, Nlat
      if (INT(Reg(ilon,ilat))==ireg) then
         if (X(ilon,ilat,it)/=Miss) then
            tmp=tmp+X(ilon,ilat,it); icount=icount+1
         endif
      endif
   Enddo
Enddo
if (icount>0) tmp=tmp/real(icount); Area_Mean = tmp
End Function Area_Mean
!--------------------------------------------------------------------------------
Subroutine Read_PDSI_nc (FName)
Character (Len=*) :: FName
Real, Allocatable, Dimension (:) :: tmplon, tmplat, tmpTime
Integer :: Nlatmp, Nlotmp, Ntmp
Character (Len=50) :: VarName

VarName="ori_pdsi"
Call rd_ncdf_3d (trim(FName), tmplon, tmplat, tmpTime, Nlotmp, Nlatmp, Ntmp, Dr_IDX_Tmp, trim(VarName))
!Allocate (Xlon(Nlotmp), Xlat(Nlatmp), Xtime(Ntmp))
!Xlon(:)=tmplon(:); Xlat(:)=tmplat(:); Xtime(:)=tmpTime(:); Ntime=Ntmp; Nlon=Nlotmp; Nlat=Nlatmp
Deallocate (tmplon, tmplat, tmpTime)

!VarName="ori_wplm"
!Call rd_ncdf_3d (trim(FName), tmplon, tmplat, tmpTime, Nlotmp, Nlatmp, Ntmp, WPLM, trim(VarName))
!Deallocate (tmplon, tmplat, tmpTime)
!VarName="ori_zind"
!Call rd_ncdf_3d (trim(FName), tmplon, tmplat, tmpTime, Nlotmp, Nlatmp, Ntmp, ZIND, trim(VarName))
!Deallocate (tmplon, tmplat, tmpTime)

!VarName="sel_pdsi"
!Call rd_ncdf_3d (trim(FName), tmplon, tmplat, tmpTime, Nlotmp, Nlatmp, Ntmp, scPDSI, trim(VarName))
!Deallocate (tmplon, tmplat, tmpTime)
!VarName="sel_wplm"
!Call rd_ncdf_3d (trim(FName), tmplon, tmplat, tmpTime, Nlotmp, Nlatmp, Ntmp, scWPLM, trim(VarName))
!Deallocate (tmplon, tmplat, tmpTime)
!VarName="sel_zind"
!Call rd_ncdf_3d (trim(FName), tmplon, tmplat, tmpTime, Nlotmp, Nlatmp, Ntmp, scZIND, trim(VarName))
!Deallocate (tmplon, tmplat, tmpTime)
End Subroutine Read_PDSI_nc
!--------------------------------------------------------------------------------
Subroutine Read_SPI_nc (FName)
Character (Len=*) :: FName
Real, Allocatable, Dimension (:) :: Lon, Lat, Time
Integer ::  Nlon1, Nlat1, Nt

!Call rd_ncdf_3d (trim(FName), Lon, Lat, Time, Nlon1, Nlat1, Nt, SPI01, "SPI_01")
!Xlon(:)=lon(:); Xlat(:)=lat(:); Xtime(:)=Time(:); Ntime=Nt; Nlon=Nlon1; Nlat=Nlat1
!Deallocate (Lon, Lat, Time)
!Call rd_ncdf_3d (trim(FName), Lon, Lat, Time, Nlon1, Nlat1, Nt, SPI03, "SPI_03")
!Deallocate (Lon, Lat, Time)
!Call rd_ncdf_3d (trim(FName), Lon, Lat, Time, Nlon1, Nlat1, Nt, SPI06, "SPI_06")
!Deallocate (Lon, Lat, Time)

Call rd_ncdf_3d (trim(FName), Lon, Lat, Time, Nlon1, Nlat1, Nt, Dr_IDX_Tmp, "SPI_12")

!Deallocate (Lon, Lat, Time)
!Call rd_ncdf_3d (trim(FName), Lon, Lat, Time, Nlon1, Nlat1, Nt, Precip, "Precip")
!Deallocate (Lon, Lat, Time)
End Subroutine Read_SPI_nc
!-----------------------------------------------------------------------------
Subroutine Ave_Std(X, N, Mean, Std)
Integer :: N, i, icount
Real :: X(N), Mean, Std
Mean=0.; Std=0.; icount=0
Do i=1, N
   if (X(i)/=Miss) then
      icount=icount+1; Mean=Mean+X(i); Std=Std+X(i)*X(i)
   endif
Enddo
if (icount>0) then
   Mean=Mean/real(icount); Std=SQRT(Std/real(icount)-Mean*Mean)
else
   Mean=Miss; Std=Miss
endif
End Subroutine Ave_Std
!-----------------------------------------------------------------------------

Include "../Libs/sub_rd_ncdf.F90"
Include "../Libs/sub_wr_ncdf_Tan.F90"
Include "All_Num_Dr_Spells.f90"  ! Tinh ca do dai tung dot

!-----------------------------------------------------------------------------
END