!-----------------------------------------------------------------------------
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
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
End Subroutine Cal_DS_DI
!-----------------------------------------------------------------------------