Program Dr_IDXs_SEACLID
  ! Calculating PDSI, SPI, SPEI 
  use netcdf

!------------------------------------------------------------------------------

  Real, parameter    :: amssng = -99., Miss=amssng
  Integer, Parameter :: NVar = 2, Nreg = 4
  Integer, Parameter :: NMonth = 12, Nmon = NMonth

  Integer, Allocatable :: Len_SPI(:)
  Real, Allocatable    :: DatBL(:,:), DatTMP(:,:), Rain(:) ! For calculating SPI, SPEI
  Real, Allocatable    :: PET(:,:)

  Real, Allocatable    :: beta_spi(:), gamm(:), pzero(:), SPI(:,:,:), SPE(:,:,:)
  Real, Allocatable, Dimension (:,:) :: mon_ori_pdsi, mon_ori_wplm, mon_ori_zind
  Real, Allocatable, Dimension (:,:) :: mon_sel_pdsi, mon_sel_wplm, mon_sel_zind

  Real, Allocatable :: Dat_Mon(:,:,:,:), DatNor(:,:,:)   ! For calculating PDSI
  Real, Allocatable :: AWC(:), xlon(:), xlat(:), Mask(:,:), xlon1(:), xlat1(:), Reg(:,:)
  Real, Allocatable :: Mask_Land(:), Lon(:), Lat(:)  ! (Nlon1*Nlat1) - All gridpoints
  !Real, Allocatable, Dimension (:,:) :: Prec, Temp
  Real, allocatable, dimension (:,:,:) :: Temp_Inp, Rain_Inp, Pet_Inp

  ! Notes:
  ! 1) Dat_Mon(ivar,ist,iyr,imw), where ivar=1:Temp (F), ivar=2:Rainfall (inch); ist=1..Nst; iyr=YBas1..YBas2; imw=1..Nmon (or Nweek) 
  !    Dat_Mon used to calculating PDSI
  ! 2) Temp_Inp (Temp in C deg), Rain_Inp (Rainfall in mm/month) (1:Nlon1,1:Nlat1,1:Ntime), Ntime=(YBas2-YBas1+1)*Nmon
  !    Temp_Inp and Rain_Inp to be used to write out as input data
  ! 3) Prec (mm/month), Temp (C) (1:Ntime,1:Nst)
  ! 4) DatBL (mm/month), DatTMP (C) (1:Ntime,1:Nst): Used to calculating SPI, SPEI

  Real, Allocatable :: T2m(:,:,:), Pr(:,:,:)

!------------------------------------------------------------------------------
  Character (Len=255) :: Fname, R_Fname, T_Fname, BL_Dir, Out_Dr_Fname, Mask_File, Out_Dir
  Character (Len=3)   :: TVar, RVar
  Character (Len=50)  :: Model_nm
  Logical :: status, Empty

  Integer :: YBas1, YBas2, NYear, maxyrs
  Integer :: NLen, Nst, Nlon, Nlat, IO, Ntime, Nt, it, Landpoint
  Integer :: imod, ist, iyr, imw, ivar, isys, system, ilon, ilat

!------------------------------------------------------------------------------
! For writing to NC file:
  Real, allocatable, Dimension (:, :, :) :: X1, X2, X3, X4, X5, X6
  Real, allocatable, Dimension (:) :: Xt
!------------------------------------------------------------------------------

  Namelist /File_Dir_Names/ R_Fname,T_Fname,Mask_File,Out_Dr_Fname,TVar, RVar 
  Namelist /Data_Info/      YBas1, YBas2, NLen
  Namelist /Dr_Timescale/   Len_SPI
!------------------------------------------------------------------------------

  Open (1, file="Namelist_GPCC_CRU_Reg.txt", Status="Old")     ! Tinh voi APH 1961-2015

  Read (1, File_Dir_Names)
  Read (1, Data_Info)
  Allocate (Len_SPI(NLen))
  Read (1, Dr_Timescale)
  
  Call Read_kmean    ! Return Reg(Nlon,Nlat), xlon(Nlon), xlat(Nlat) 
  Nyear = (YBas2 - YBas1 + 1); Ntime = Nyear*NMonth; maxyrs = NYear
  Nst = Nreg + 1  ! 4 sub-reg + 1 entire domain

  print*,"Nst,Ntime,NLen=",Nst,Ntime,NLen
  Allocate(Temp_Inp(1:Nlon,1:Nlat,1:Ntime));
  Allocate(Rain_Inp(1:Nlon,1:Nlat,1:Ntime));
  Allocate(Pet_Inp(1:Nlon,1:Nlat,1:Ntime))      ! Boc hoi 

     Allocate (Rain(Ntime),DatBL(Ntime,Nst),DatTMP(Ntime,Nst))
     Allocate (PET(Ntime,Nst))
     Allocate (beta_spi(Nmonth), gamm(Nmonth), pzero(Nmonth))
     Allocate (SPI(Nst,Ntime,NLen), SPE(Nst,Ntime,NLen))

     !R_Fname="/home/tanpv/Data/T_R_APH_CRU_GPCC/T_P_Merge_61-15.nc"
     !T_Fname=R_Fname

     Call Read_Pr_T2m_Merge_Data
     print*,"Reading Data: Done: Nlat,Nlon,Nst=",Nlat,Nlon,Nst

     ! Calculating SPI:
     !------------------------------------------------------------------------------
     Do it=1,Ntime; print"(i4,20F7.1)",it, DatTMP(it,:),DatBL(it,:),PET(it,:); Enddo;pause
     SPI(:,:,:) = amssng; SPE(:,:,:) = amssng
     Call Cal_SPEI_Reg              ! <<<<<<<<  Should check again if SPE is SPI or SPEI  ????
     SPE(:,:,:) = SPI(:,:,:)

     ! Calculating SPI:
     !------------------------------------------------------------------------------
     SPI(:,:,:) = amssng
     Call Cal_SPI

     print*,"Calculating SPI and SPEI done !"
     !Fname=trim(Out_Dir)//"BL_"//trim(Model_nm)//"_SEA_Reg_61-15.nc"
     Call Write_All_Dr_IDX_to_NetCDF (trim(Out_Dr_Fname))
CONTAINS
!------------------------------------------------------------------------------
Subroutine Read_kmean 
   Integer :: ist, ilon, ilat, nlo, nla
   Real, Allocatable :: xtmp(:,:),lontmp(:),lattmp(:)
   Character (Len=255) :: Fname

     Fname="/home/tanpv/Progs/2021/NC/kmeans_4_modified.nc"
     Call rd_ncdf_2d (Fname,xlon,xlat,Nlon,Nlat,Reg, "cluster")  ! Sub-Region
     Where (Reg==-9999.0) Reg = Miss
End Subroutine Read_kmean 
   !------------------------------------------------------------------------------
   Subroutine Fill_up (X,Ntime)
      Real :: X(Ntime), TB(Nmon)
      Integer :: ic, it, imon, imon1, Nleng
      Real :: tmp

      if (mod(Ntime,Nmon)/=0) print*,"Something wrong !"
      Nleng=Int(0.8*Ntime)   ! 80% do dai chuoi co SL
      Do imon=1,Nmon
         tmp=0; ic=0
         Do it=1,Ntime
            imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon
            if (imon1==imon) then
               if (X(it)/=Miss) then
                  ic=ic+1; tmp=tmp+X(it)
               endif
            endif
         Enddo
         if (ic>=Nleng) then; TB(imon)=tmp/ic; else; TB(imon)=Miss; endif
      Enddo
      Do it=1,Ntime
         imon=mod(it,Nmon); if (imon==0) imon=Nmon
         if (TB(imon)/=Miss.AND.X(it)==Miss) X(it)=TB(imon)
      Enddo
   End Subroutine Fill_up
   !------------------------------------------------------------------------------
   Logical Function FullDara(X,N,Miss)
      Integer :: N, i
      Real :: X(N), Miss
      Logical :: Full

      Full=.True.
      Do i=1,N
         if (X(i)==Miss) then; Full=.False.; exit; endif
      Enddo
      FullDara=Full
   End Function FullDara
   !------------------------------------------------------------------------------
   Subroutine Read_GPCC_CRU
      Integer, Parameter :: Yr1=1951, Yr2=2019
      Real, Parameter :: Miss_gpcc=-99999.99, Miss_cru=9.96921e+36
      Real, Allocatable :: xtime(:), T2m(:,:,:), Pr(:,:,:), x1d(:), y1d(:)
      Real, Allocatable :: xlon(:), xlat(:)
      Integer :: Nt, Nt1, Nt2, ilon, ilat, it, it1, it2
      ! xlon(:), xlat(:), Mask(:,:)

      Call rd_ncdf_3d (trim(R_Fname), xlon, xlat, xtime, Nlon, Nlat, Nt1, Pr, trim(RVar))
      Deallocate (xlon, xlat, xtime)
      Call rd_ncdf_3d (trim(T_Fname), xlon, xlat, xtime, Nlon, Nlat, Nt2, T2m, trim(TVar))
      Deallocate (xlon, xlat)
      Call rd_ncdf_2d (trim(Mask_File),xlon,xlat,nlon,nlat,Mask,"land")

      Where (Pr==Miss_gpcc) Pr=Miss
      Where (T2m==Miss_cru) T2m=Miss

      it1=(YBas1-Yr1)*NMonth+1  ! 1960-01
      it2=it1+Ntime-1           ! 2015-12
      print*,"Nlon,Nlat,Ntime,Nt1,Nt2,it1,it2=",Nlon,Nlat,Ntime,Nt1,Nt2,it1,it2

      !Allocate (Temp_Inp(Nlon,Nlat,Ntime),Rain_Inp(Nlon,Nlat,Ntime))
      Allocate (x1d(Ntime),y1d(Ntime))
      Temp_Inp=Miss; Rain_Inp=Miss
      Do ilon=1,Nlon
         Do ilat=1,Nlat
            if (Mask(ilon,ilat)/=Miss) then
               x1d(1:Ntime)=Pr(ilon,ilat,it1:it2)
               Call Fill_up (x1d,Ntime)
               Rain_Inp(ilon,ilat,1:Ntime)=x1d(1:Ntime)

               x1d(1:Ntime)=T2m(ilon,ilat,it1:it2)
               Call Fill_up (x1d,Ntime)
               Temp_Inp(ilon,ilat,1:Ntime)=x1d(1:Ntime)
            endif
         Enddo
      Enddo
      Do ilon=1,Nlon
         Do ilat=1,Nlat
            if (Mask(ilon,ilat)/=Miss) then
               x1d(1:Ntime)=Rain_Inp(ilon,ilat,1:Ntime)
               y1d(1:Ntime)=Temp_Inp(ilon,ilat,1:Ntime)
               if (.not.FullDara(x1d,Ntime,Miss).AND..not.FullDara(y1d,Ntime,Miss)) then
                  Mask(ilon,ilat)=Miss
                  Rain_Inp(ilon,ilat,1:Ntime)=Miss
                  Temp_Inp(ilon,ilat,1:Ntime)=Miss
               endif
            endif
         Enddo
      Enddo
      !Nst=0
      !Do ilat=1,Nlat
      !   Do ilon=1,Nlon
      !      if (Mask(ilon,ilat)/=Miss) Nst=Nst+1
      !   Enddo
      !Enddo
   End Subroutine Read_GPCC_CRU
!------------------------------------------------------------------------------
Subroutine Read_Pr_T2m_Merge_Data
  ! Read data from file in which data are merged from APH, CRU and GPCC
  ! Data from 1961-2005
  !Real, Allocatable :: XTime(:)
  !Real, Allocatable :: xlon(:), xlat(:), ET(:), Temp(:)
  Real, Allocatable :: ET(:), Temp(:)
  Integer :: Nc1, Nc2, Nc3, ireg, Nlon1, Nlat1, Nt0, it1, it2, it, it0,ilon,ilat
  Real :: tmp1, tmp2, tmp3

  !Call rd_ncdf_3d (trim(R_Fname), xlon, xlat, XTime, Nlon1, Nlat1, Nt, Pr,  trim(RVar))
  !Deallocate (xlon, xlat, XTime)
  !Call rd_ncdf_3d (trim(T_Fname), xlon, xlat, XTime, Nlon1, Nlat1, Nt, T2m, trim(TVar))

  Call Read_GPCC_CRU

  Allocate (ET(0:Ntime-1), temp(0:Ntime-1))
  Do ilat=1, Nlat     ! Only grids over SEA (>=12 deg N) 
     Do ilon=1, Nlon
        if (Reg(ilon,ilat)/=Miss) then
           !Temp_Inp(ilon,ilat,1:Ntime)=T2m(ilon,ilat,1:Ntime)
           !Rain_Inp(ilon,ilat,1:Ntime)=Pr (ilon,ilat,1:Ntime)
           temp(0:Ntime-1)=Temp_Inp(ilon,ilat,1:Ntime)
           Call PET_thornthwaite(temp,Ntime,xlat(ilat),ET)
           Pet_Inp(ilon,ilat,1:Ntime) = ET(0:Ntime-1) 
        else
           Temp_Inp(ilon,ilat,1:Ntime)=Miss
           Rain_Inp(ilon,ilat,1:Ntime)=Miss
           Pet_Inp(ilon,ilat,1:Ntime) =Miss
        endif
     Enddo
  Enddo
  ! Areal mean:
  ! DatBL(Ntime,Nst),DatTMP(Ntime,Nst)
  DatBL = Miss; DatTMP = Miss
  Do it=1,Ntime
     ! Sub-region
     Do ireg=1,Nreg
        Nc1=0; tmp1=0; Nc2=0; tmp2=0; Nc3=0; tmp3=0
        Do ilon=1,Nlon
           Do ilat=1,Nlat
              if (Reg(ilon,ilat)==ireg) then
                 if (Temp_Inp(ilon,ilat,it)/=Miss) then
                    Nc1=Nc1+1; tmp1=tmp1+Temp_Inp(ilon,ilat,it)
                 endif
                 if (Rain_Inp(ilon,ilat,it)/=Miss) then
                    Nc2=Nc2+1; tmp2=tmp2+Rain_Inp(ilon,ilat,it)
                 endif
                 if (Pet_Inp(ilon,ilat,it)/=Miss) then
                    Nc3=Nc3+1; tmp3=tmp3+Pet_Inp(ilon,ilat,it)
                 endif
              endif
           Enddo
        Enddo
        if (Nc1>0) then; DatTMP(it,ireg)=tmp1/Nc1; else; DatTMP(it,ireg)=Miss; endif
        if (Nc2>0) then; DatBL (it,ireg)=tmp2/Nc2; else; DatBL (it,ireg)=Miss; endif
        if (Nc3>0) then; PET (it,ireg)=tmp3/Nc3; else; PET (it,ireg)=Miss; endif
     Enddo
     ! Entire region
     Nc1=0; tmp1=0; Nc2=0; tmp2=0; Nc3=0; tmp3=0
     Do ilon=1,Nlon
        Do ilat=1,Nlat
           if (Reg(ilon,ilat)>0) then
              if (Temp_Inp(ilon,ilat,it)/=Miss) then
                 Nc1=Nc1+1; tmp1=tmp1+Temp_Inp(ilon,ilat,it)
              endif
              if (Rain_Inp(ilon,ilat,it)/=Miss) then
                 Nc2=Nc2+1; tmp2=tmp2+Rain_Inp(ilon,ilat,it)
              endif
              if (Pet_Inp(ilon,ilat,it)/=Miss) then
                 Nc3=Nc3+1; tmp3=tmp3+Pet_Inp(ilon,ilat,it)
              endif
           endif
        Enddo
     Enddo
     if (Nc1>0) then; DatTMP(it,Nreg+1)=tmp1/Nc1; else; DatTMP(it,Nreg+1)=Miss; endif
     if (Nc2>0) then; DatBL (it,Nreg+1)=tmp2/Nc2; else; DatBL (it,Nreg+1)=Miss; endif
     if (Nc3>0) then; PET (it,Nreg+1)=tmp3/Nc3; else; PET (it,Nreg+1)=Miss; endif
     !print"(I4,20F8.1)",it,(DatBL(it,i),DatTMP(it,i),DatBL(it,i)-PET(it,i),i=1,Nreg+1)
  Enddo
  !pause
End Subroutine Read_Pr_T2m_Merge_Data
!------------------------------------------------------------------------------
Subroutine Write_All_Dr_IDX_to_NetCDF(Fname) !, Nlon, Nlat, Ntime)
   Integer :: imon, imon1, ist, jlat, ilon, it, ilen, ix
   Character (Len=50) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, t_nm
   Character (Len=*) :: Fname
   Real, Allocatable :: xreg(:), xtime(:), xlen(:), X1(:,:,:), X2(:,:,:), X3(:,:,:), X4(:,:,:),X5(:,:,:)
   Real, Allocatable :: X6(:,:,:), X7(:,:,:),X8(:,:,:)
   Real, Allocatable :: Pr_ano(:,:), T_ano(:,:), Pet_ano(:,:)
   Real :: tbr(Nmon), tbt(Nmon), tbp(Nmon), tmp1,tmp2,tmp3,ntmp1,ntmp2,ntmp3 

   ! SPI(Nst,Ntime,NLen), SPE(Nst,Ntime,NLen))
   ! Allocate (Rain(Ntime),DatBL(Ntime,Nst),DatTMP(Ntime,Nst))
   ! Allocate (PET(Ntime,Nst))

   ! Tinh Anomaly cua Rain, Temp, Pet:
   Allocate (Pr_ano(Ntime,Nst),T_ano(Ntime,Nst),Pet_ano(Ntime,Nst))
   Pr_ano=Miss; T_ano=Miss; Pet_ano=Miss
      ! Monthly climatology:
   Do ist=1,Nst
      Do imon=1,Nmon  ! Monthly climatology
         tmp1=0; tmp2=0; tmp3=0; ntmp1=0; ntmp2=0; ntmp3=0
         Do it=1,Ntime
            imon1=Mod(it,Nmon); if (imon1==0) imon1=Nmon
            if (DatBL(it,ist)/=Miss) then
               ntmp1=ntmp1+1; tmp1=tmp1+DatBL(it,ist)
            endif
            if (DatTMP(it,ist)/=Miss) then
               ntmp2=ntmp2+1; tmp2=tmp2+DatTMP(it,ist)
            endif
            if (PET(it,ist)/=Miss) then
               ntmp3=ntmp3+1; tmp3=tmp3+PET(it,ist)
            endif
         Enddo
         if (ntmp1/=0) then; tbr(imon)=tmp1/ntmp1; else; tbr(imon)=Miss; endif
         if (ntmp2/=0) then; tbt(imon)=tmp2/ntmp2; else; tbt(imon)=Miss; endif
         if (ntmp3/=0) then; tbp(imon)=tmp3/ntmp3; else; tbp(imon)=Miss; endif
      Enddo
      Do imon=1,Nmon  ! Anomaly 
         Do it=1,Ntime
            imon1=Mod(it,Nmon); if (imon1==0) imon1=Nmon
            if (tbr(imon1)/=Miss) Pr_ano (it,ist)=DatBL (it,ist)-tbr(imon1)
            if (tbt(imon1)/=Miss) T_ano  (it,ist)=DatTMP(it,ist)-tbt(imon1)
            if (tbp(imon1)/=Miss) Pet_ano(it,ist)=PET   (it,ist)-tbp(imon1)
         Enddo
      Enddo
   Enddo
   ! Writing:
   Allocate (xreg(Nst),xtime(Ntime),xlen(Nlen))
   Allocate (X1(Nst,Nlen,Ntime),X2(Nst,Nlen,Ntime))
   Allocate (X3(Nst,Nlen,Ntime),X4(Nst,Nlen,Ntime),X5(Nst,Nlen,Ntime))
   Allocate (X6(Nst,Nlen,Ntime),X7(Nst,Nlen,Ntime),X8(Nst,Nlen,Ntime))

   Do ix=1,Nst; xreg(ix)=ix; Enddo
   Do it=1,Ntime; xtime(it)=it; Enddo
   Do ix=1,Nlen; xlen(ix)=ix; Enddo

   Do ix=1,Nst
      Do ilen=1,Nlen
         X1(ix,ilen,1:Ntime)=SPI(ix,1:Ntime,ilen)   ! spi
         X2(ix,ilen,1:Ntime)=SPE(ix,1:Ntime,ilen)   ! spei

         X6(ix,ilen,1:Ntime)=DatBL(1:Ntime,ix)      ! Rainfall
         X7(ix,ilen,1:Ntime)=DatTMP(1:Ntime,ix)     ! T2m
         X8(ix,ilen,1:Ntime)=PET(1:Ntime,ix)        ! PET

         X3(ix,ilen,1:Ntime)=Pr_ano (1:Ntime,ix)      ! Rainfall Ano
         X4(ix,ilen,1:Ntime)=T_ano  (1:Ntime,ix)      ! T2m Ano
         X5(ix,ilen,1:Ntime)=Pet_ano(1:Ntime,ix)     ! PET Ano
      Enddo
   Enddo

   Write(t_unit,"(A,I4.4,A)") "months since ",YBas1,"-01-01"
   Var_ln="spi-12"; var_unit="-"; lon_nm="reg"; lat_nm="timescale"; t_nm="month"
   
   Call wr_ncdf_3d (Fname,xreg,xlen,xtime,Nst,Nlen,Ntime,X1,"spi",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
   Var_ln="spei-12";
   Call wr_ncdf_3d_apd (Fname,X2,"spe12",Nst,Nlen,Ntime, Var_ln, var_unit,lon_nm,lat_nm,t_nm,Miss)
   !Var_ln="Monthly rainfall averaged over the sub-reg"; var_unit="mm/mon"
   Var_ln="Monthly Anomaly rainfall over the sub-reg"; var_unit="mm/mon"
   Call wr_ncdf_3d_apd (Fname,X3,"pra",Nst,Nlen,Ntime, Var_ln, var_unit,lon_nm,lat_nm,t_nm,Miss)
   !Var_ln="Monthly 2m-temperature averaged over the sub-reg"; var_unit="oC"
   Var_ln="Monthly Anomaly 2m-temperature over the sub-reg"; var_unit="oC"
   Call wr_ncdf_3d_apd (Fname,X4,"t2ma",Nst,Nlen,Ntime, Var_ln, var_unit,lon_nm,lat_nm,t_nm,Miss)
   !Var_ln="Monthly potential evaporation calculated for each sub-reg"; var_unit="mm/mon"
   Var_ln="Monthly Anomaly potential evaporation for each sub-reg"; var_unit="mm/mon"
   Call wr_ncdf_3d_apd (Fname,X5,"peta",Nst,Nlen,Ntime, Var_ln, var_unit,lon_nm,lat_nm,t_nm,Miss)

   Var_ln="Monthly rainfall averaged over the sub-reg"; var_unit="mm/mon"
   Call wr_ncdf_3d_apd (Fname,X6,"pr",Nst,Nlen,Ntime, Var_ln, var_unit,lon_nm,lat_nm,t_nm,Miss)
   Var_ln="Monthly 2m-temperature averaged over the sub-reg"; var_unit="oC"
   Call wr_ncdf_3d_apd (Fname,X7,"t2m",Nst,Nlen,Ntime, Var_ln, var_unit,lon_nm,lat_nm,t_nm,Miss)
   Var_ln="Monthly potential evaporation calculated for each sub-reg"; var_unit="mm/mon"
   Call wr_ncdf_3d_apd (Fname,X8,"pet",Nst,Nlen,Ntime, Var_ln, var_unit,lon_nm,lat_nm,t_nm,Miss)
End Subroutine Write_All_Dr_IDX_to_NetCDF
!------------------------------------------------------------------------------
Logical Function Empty_Stn_Arr1 (X,N)
   ! Tra ve .True. neu tram ist khuyet SL
   Integer :: N
   Real :: X(N)
   Integer :: i,iyr,imon, k1, k2
   Logical :: Empty
     Empty=.false.
     Do i=1,N
        if (X(i)==Miss) then
           Empty=.true.; exit
        endif
     Enddo
     Empty_Stn_Arr1 = Empty   !    ! T hoac P khuyet SL
End Function Empty_Stn_Arr1
!------------------------------------------------------------------------------
Logical Function Empty_Stn (DatBL, ist)
   ! Tra ve .True. neu tram ist khuyet SL
   Real, Dimension (NVar, Nst, YBas1:YBas2, NMonth) :: DatBL
   Integer :: N, ist,iyr,imon, k1, k2
   Logical :: L1, L2

   L1=.False.; N=0
   Do iyr=YBas1, YBas2
      Do imon=1,NMonth
         if (DatBL(1,ist,iyr,imon)<0.OR.DatBL(2,ist,iyr,imon)<0) then
            N = N+1
            if (N>10) then ; L1=.True.; Exit; endif 
         endif
      Enddo
   Enddo
   Empty_Stn = L1 !.OR.L2    ! T hoac P khuyet SL
End Function Empty_Stn
!------------------------------------------------------------------------------
  Integer Function DayOfMon(Mon, Year)
  Integer :: Mon, Year
  Select Case (Mon)
     Case (1,3,5,7,8,10,12)
        DayOfMon = 31
     Case (4,6,9,11)
        DayOfMon = 30
     Case (2)
        if (mod(Year,4)==0.AND.mod(Year,100)/=0) then
           DayOfMon = 29
        else
           DayOfMon = 28
        endif
  End Select
  End Function DayOfMon
!--------------------------------------------------------------
Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"

Include "/home/tanpv/Mylibs/Sub_SPI.f90"  ! Chuong trinh tinh chi so han SPI
Include "/home/tanpv/Mylibs/Sub_SPEI_Reg.f90" ! Chuong trinh tinh chi so han SPEI
Include "/home/tanpv/Mylibs/PET_thornthwaite.f90"
Include "/home/tanpv/Mylibs/Percentile.f90"

!--------------------------------------------------------------------
END ! Program