! 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, N_Drou_IDX
   Dr_IDX=trim(Drou_IDX(iidx))
   Do ifile=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
!-----------------------------------------------------------------------------
Include "SUB/Area_Mean.f90"
Include "SUB/Ave_Std.f90"
Include "SUB/Cal_DS_DI.f90"
Include "SUB/Cal_Freq_All_Timeseries.f90"
Include "SUB/Cal_Freq_Monthly.f90"
Include "SUB/Cal_Spells.f90"
Include "SUB/Define_Drought_Events.f90"
Include "SUB/Drought_Index_Areal_Mean.f90"
Include "SUB/Geographic_Extent_1.f90"
!Include "SUB/Geographic_Extent.f90"
Include "SUB/Read_Data.f90"
Include "SUB/Reg_Drought_Event.f90"
Include "SUB/Write_Freq_Duration_To_NC.f90"
Include "SUB/Write_Severity_Intensity_To_NC.f90"
!Include "SUB/Write_Spell_Area_To_TXT.f90"
Include "SUB/Write_Spell_GeoExtent_To_NC_1.f90"
Include "SUB/Write_Spell_To_NC.f90"

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