Program Dr_IDXs_SEACLID
  ! Calculating PDSI, SPI, SPEI 
  use netcdf

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

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

  Integer, Allocatable :: Len_SPI(:)
  Real, Allocatable    :: DatBL(:,:), DatTMP(:,:), Rain(:)

  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, Dimension (:,:) :: Rain_Mon

  Real, Allocatable :: Dat_Mon(:,:,:,:), DatNor(:,:,:)
  Real, Allocatable :: AWC(:), xlon(:), xlat(:), Mask(:,:), Mask_Land(:), Lon(:), Lat(:)
  Real, Allocatable, Dimension (:,:) :: Prec, Temp
  Real, allocatable, dimension (:,:,:) :: Temp_Inp, Rain_Inp


!------------------------------------------------------------------------------
  Character (Len=255) :: Fname, BL_Dir, PDSI_Dir, PDSI_WDir, R_Fname, T_Fname, MaskFname, AWC_Fname, Out_Dr_Fname
  Character (Len=3)   :: TVar, RVar
  Character (Len=30) :: FMT
  Logical :: status, Empty

  Integer :: YBas1, YBas2, NYear, maxyrs
  Integer :: NLen, Nst, Nlon, Nlat, IO, Ntime, it
  Integer :: 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/ PDSI_Dir, PDSI_WDir, R_Fname, T_Fname, MaskFname, AWC_Fname, Out_Dr_Fname, &
                            TVar, RVar
  Namelist /Data_Info/      YBas1, YBas2, NLen
  Namelist /Dr_Timescale/   Len_SPI
!------------------------------------------------------------------------------

  Open (1, file="Namelist_APH.txt", Status="Old")
  Read (1, File_Dir_Names)
  Read (1, Data_Info)
  Allocate (Len_SPI(NLen))
  Read (1, Dr_Timescale)

  Nyear = (YBas2 - YBas1 + 1); Ntime = Nyear*NMonth; maxyrs = NYear
  Call Read_APH_netcdf 
  print*,"Reading Data: Done !"

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

  Allocate (DatNor (NVar,Nst,Nmonth), Rain_Mon(Nst,Ntime))
  Allocate (AWC(Nst)) !, xlon(Nst), xlat(Nst))
  Allocate (mon_ori_pdsi(Nst,Ntime), mon_ori_wplm(Nst,Ntime), mon_ori_zind(Nst,Ntime))
  Allocate (mon_sel_pdsi(Nst,Ntime), mon_sel_wplm(Nst,Ntime), mon_sel_zind(Nst,Ntime))

  DatNor(:,:,:) = 0.; Rain_Mon(:,:)=Miss !-99

  Do ist=1, Nst ! For calculating SPI and SPEI
     DatBL(1:Ntime,ist)  = Prec(ist,1:Ntime);
     DatTMP(1:Ntime,ist) = Temp(ist,1:Ntime);
  Enddo

  ! Calculating PDSI:
  !------------------------------------------------------------------------------
  Do ist=1, Nst
     it=0
     Do iyr=YBas1, YBas2
        Do imw=1, NMonth
           it=it+1; Rain_Mon(ist,it) = Dat_Mon(2,ist,iyr,imw)
        Enddo
     Enddo
  Enddo

  Call putAWC(AWC, lon, lat, Nst) !; Where (AWC<0) AWC=6.7 ! Assign to gridboxes located along coast
  isys = system ("rm -Rf "//trim(PDSI_Dir))
  isys = system ("mkdir -p "//trim(PDSI_Dir))
  mon_ori_pdsi(:,:)=Miss; mon_ori_wplm(:,:)=Miss; mon_ori_zind(:,:)=Miss;
  mon_sel_pdsi(:,:)=Miss; mon_sel_wplm(:,:)=Miss; mon_sel_zind(:,:)=Miss;

  Allocate (Mask_Land(Nst))
  ist=0
  Do ilat=1, Nlat; Do ilon=1, Nlon; ist=ist+1; Mask_Land(ist)=Mask(ilon,ilat); Enddo; Enddo
  Do ist=1, Nst
  !   print*,AWC(ist), lat(ist)
     if (Mask_Land(ist)>0) then
        if (AWC(ist)<0) AWC(ist)=6.7
        Print*, "Calculation for station: ", ist  
  !      !Call Write_to_PDSI ("monthly",  Dat_Mon,  NMonth, ist)
        Call Write_to_PDSI ("monthly",NMonth, ist)
        isys = system ("./pdsi.exe ")
  !      ! PostProcessing for station
        Call Read_Table (trim(PDSI_Dir)//"/monthly/original/"//"PDSI.tbl", NMonth, mon_ori_pdsi(ist,:))
        Call Read_Table (trim(PDSI_Dir)//"/monthly/original/"//"WPLM.tbl", NMonth, mon_ori_wplm(ist,:))
        Call Read_Table (trim(PDSI_Dir)//"/monthly/original/"//"ZIND.tbl", NMonth, mon_ori_zind(ist,:))
        Call Read_Table (trim(PDSI_Dir)//"/monthly/self_cal/"//"PDSI.tbl", NMonth, mon_sel_pdsi(ist,:))
        Call Read_Table (trim(PDSI_Dir)//"/monthly/self_cal/"//"WPLM.tbl", NMonth, mon_sel_wplm(ist,:))
        Call Read_Table (trim(PDSI_Dir)//"/monthly/self_cal/"//"ZIND.tbl", NMonth, mon_sel_zind(ist,:))
     endif
  Enddo
  print*,"Calculating PDSI: Done !"

  ! Calculating SPI:
  !------------------------------------------------------------------------------
  SPI(:,:,:) = amssng; SPE(:,:,:) = amssng
  Call Cal_SPEI              ! <<<<<<<<  Should check again if SPE is SPI or SPEI  ????
  SPE(:,:,:) = SPI(:,:,:)

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

  Call Write_All_Dr_IDX_to_NetCDF (trim(Out_Dr_Fname))

CONTAINS
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
Subroutine Read_Table (Fname, Nmon_week, Outp) 
   Integer :: Nmon_week, it, iyr, imw, year, isys
   Real :: wk(Nmon_week), Outp(Nmon_week*NYear)
   Character (Len=*) :: Fname
   Logical :: exists

   Inquire (File=trim(Fname),exist=exists)
   if (.Not.exists) return
   Open (1, File=trim(Fname), status="Old")
   it = 0
   !print*,trim(Fname), "<<<<<<================"
   Do iyr=YBas1, YBas2
      read (1, *) year, (wk(imw), imw=1, Nmon_week)
      !   print*, year, (wk(imw), imw=1, Nmon_week)
      Do imw=1, Nmon_week
         it = it+1; Outp(it) = wk(imw)
      Enddo
   Enddo
   close (1)
   isys = system ("rm -f "//trim(Fname))
End Subroutine Read_Table
!------------------------------------------------------------------------------
!Subroutine Write_to_PDSI (Dat_Dir, DatBL, Nmon_week, ist)

Subroutine Write_to_PDSI (Dat_Dir,Nmon_week, ist)

   Integer :: ivar, ist, Nmon_week, iyr, imw, isys, system, ic
   Real :: Tmp
   Character (Len=*) :: Dat_Dir   ! Dat_Dir="monthly"
   Character (Len=155) :: Fname
   !Real, Dimension (NVar, Nst, YBas1:YBas2, NMonth) :: DatBL

   ! Write data of one station to calculate PDSI
   Fname = PDSI_WDir !"/home/tanpv/SEACLID_New/PDSI/my_dat"
   !print*,ist,trim(Fname)
   isys = system ('mkdir -p '//trim(Fname))
   !print*,"mkdir OK"
   !print*,"File=",trim(Fname)//"/"//"parameter"
   !Open (1, file=trim(Fname)//"/"//"parameter", status="Unknown")
   !print*,"Open file OK", AWC(ist), xlat(ist)
   !write (1,"(2F20.5)") AWC(ist), xlat(ist)
   !close (1) 
   !print*,"AWC,Xlat= >>>>>>" !,AWC(ist), xlat(ist); 
   Do ivar=1, NVar
      if (ivar == 1) then
         Open (3, file=trim(Fname)//"/"//trim(Dat_Dir)//"_T", status="Unknown")
         Open (4, file=trim(Fname)//"/"//Dat_Dir(1:3)//"_T_normal",status="Unknown")
      else
         Open (3, file=trim(Fname)//"/"//trim(Dat_Dir)//"_P", status="Unknown")
         Open (4, file=trim(Fname)//"/"//Dat_Dir(1:3)//"_P_normal",status="Unknown")
      endif
      Do imw=1, Nmon_week   ! Tinh trung binh khi hau
         Tmp = 0; ic=0
         Do iyr=YBas1, YBas2
            !Tmp = Tmp + DatBL(ivar,ist,iyr,imw)
            if (Dat_Mon(ivar,ist,iyr,imw)/=Miss) then
               Tmp = Tmp + Dat_Mon(ivar,ist,iyr,imw)
               ic=ic+1                                ! So nam co SL hop le
            endif
            !print*,Dat_Mon(ivar,ist,iyr,imw),Tmp
         Enddo
         DatNor(ivar,ist,imw) = Tmp/Real(ic)
      Enddo
      write (4,"(52F8.3)") (DatNor(ivar,ist,imw), imw=1, Nmon_week)
      Do iyr=YBas1, YBas2
         !write (3,"(I4,52F8.3)") iyr,(DatBL(ivar,ist,iyr,imw), imw=1, Nmon_week)
         Do imw=1, Nmon_week
            if (Dat_Mon(ivar,ist,iyr,imw)==Miss) Dat_Mon(ivar,ist,iyr,imw)=DatNor(ivar,ist,imw)   ! Thay bang TB nhieu nam
         Enddo
         write (3,"(I4,52F8.3)") iyr,(Dat_Mon(ivar,ist,iyr,imw), imw=1, Nmon_week) 
!         write (*,"(I4,52F8.3)") iyr,(DatBL(ivar,ist,iyr,imw), imw=1, Nmon_week)
      Enddo
      Close (3); Close (4)
   Enddo
   Fname=trim(PDSI_WDir)//"/"//"parameter"
   Open (Unit=12, File=trim(Fname), Status="Unknown") !,FORM="FORMATTED")   
   write (12,"(2F20.5)") AWC(ist), lat(ist)
   Close (12)
End Subroutine Write_to_PDSI
!------------------------------------------------------------------------------
Subroutine Read_APH_netcdf !(T_Fname,R_Fname,TVar,RVar)
!  Character (Len=*) :: R_Fname,T_Fname,TVar,RVar
  Real, Allocatable :: XTime(:), T2m(:,:,:), Pr(:,:,:) !, xlon1(:), xlat1(:) !, VN(:,:)
  Integer :: nxlon, nxlat, ilon, jlat
  Integer :: Nlon1, Nlat1, Nt, i, j, ist, it, iyr,imw, ict,icr, imon1
  Integer, parameter :: NMon_Skip = 72    ! Bo qua 6 Yrs x 12 mon = 72 mons tu 1970-1975 (Lay tu 1976 - 2005)
  Real :: TMean(Nmonth),RMean(Nmonth)

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

! Luu Temperature and Precipitation for output:
  if (.Not.Allocated(Temp_Inp)) Allocate(Temp_Inp(1:Nlon1,1:Nlat1,1:Ntime)); 
  if (.Not.Allocated(Rain_Inp)) Allocate(Rain_Inp(1:Nlon1,1:Nlat1,1:Ntime));
  if (.Not.Allocated(Mask))     Allocate(Mask (1:Nlon1,1:Nlat1))

  Mask(:,:) = 1
  Where (T2m(:,:,NMon_Skip)==-99.9) Mask = Miss
  Where (Pr (:,:,NMon_Skip)==-99.9) Mask = Miss
 
  Where(T2m==-99.9) T2m=Miss; Temp_Inp(:,:,1:Ntime)=T2m(:,:,NMon_Skip+1:Nt); 
  Where(Pr ==-99.9) Pr =Miss; Rain_Inp(:,:,1:Ntime)=Pr (:,:,NMon_Skip+1:Nt)
 
  ! Fill up missing data:
  Do j=1, Nlat1
     Do i=1, Nlon1
        if (Mask(i,j) == 1) then
           Do imon=1, Nmonth
              TMean(imon)=0; RMean(imon)=0; ict=0; icr=0
              Do it=1, Ntime
                 imon1=mod(it,Nmonth); if (imon1==0) imon1=Nmonth
                 if (imon1==imon) then
                    if (Temp_Inp(i,j,it)/=Miss) then 
                       TMean(imon)=TMean(imon)+Temp_Inp(i,j,it); ict=ict+1
                    endif
                    if (Rain_Inp(i,j,it)/=Miss) then
                       RMean(imon)=RMean(imon)+Rain_Inp(i,j,it); icr=icr+1
                    endif
                 endif
              Enddo
              TMean(imon)=TMean(imon)/Real(ict)
              RMean(imon)=RMean(imon)/Real(icr)
           Enddo
           Do it=1, Ntime
              imon=mod(it,Nmonth); if (imon==0) imon=Nmonth
              if (Temp_Inp(i,j,it)==Miss) Temp_Inp(i,j,it)=TMean(imon)
              if (Rain_Inp(i,j,it)==Miss) Rain_Inp(i,j,it)=RMean(imon)
           ENddo
        endif
     Enddo
  Enddo

  Deallocate(T2m, Pr); Allocate (T2m(1:Nlon1,1:Nlat1,1:Ntime),Pr(1:Nlon1,1:Nlat1,1:Ntime))
  T2m = Temp_Inp; Pr = Rain_Inp   ! C and mm/day


  !print*,"Rainfall:",Pr;pause
 
  print*,"In Read_data_netcdf: Nlon, Nlat, Nt=", Nlon1, Nlat1, Ntime
  print*,lon, lat

  !Where (Pr < 0) Pr = -99; Where (T2m < 0) T2m = -99; ist = 0
  Allocate (lon(Nlon1*Nlat1), lat(Nlon1*Nlat1))
  Allocate (Dat_Mon(NVar,Nlon1*Nlat1,YBas1:YBas2, Nmonth))
  Dat_Mon(:,:,:,:) = Miss
  ist=0
  Do j=1, Nlat1
     Do i=1, Nlon1
        !if (Mask(i,j) == 1) then
           ist=ist+1; !write(Station(ist),"(A,I4.4)") "Point-",ist
           Lon(ist)=xlon(i); Lat(ist)=xlat(j);
           it=0
           Do iyr=YBas1, YBas2  ! 1, NYear
              Do imw=1, Nmonth
                 it=it+1
                 if (T2m(i,j,it)/=Miss) Dat_Mon(1,ist,iyr,imw) = T2m(i,j,it)*1.8 + 32.   ! Doi sang do F
                 if (Pr (i,j,it)/=Miss) then
                    Dat_Mon(2,ist,iyr,imw) = Pr (i,j,it)*DayOfMon(imw,iyr)/25.4          ! Doi sang Inch/month
                    !print*,DayOfMon(imw,iyr), Dat_Mon(2,ist,iyr,imw)
                    Rain_Inp(i,j,it) = Pr (i,j,it)*DayOfMon(imw,iyr)                     ! in mm/month
                 endif
              Enddo
           Enddo
        !endif
     Enddo
  Enddo

  !print*,Dat_Mon(1,:,:,:); pause
  !print*,"Rainfall:",NVar,Dat_Mon; pause 
  ! Nst = ist = Nlon*Nlat = So diem luoi trong mien

  Nst = ist; Nlon = Nlon1; Nlat = Nlat1
  Allocate (Prec(Nst,Ntime), Temp(Nst,Ntime))
  !Xlon(1:Nst)=Xlon1(1:Nst); Xlat(1:Nst)=Xlat1(1:Nst)
  ist=0
  Do j=1, Nlat !i=1, Nlon
     Do i=1, Nlon !j=1, Nlat
        ist=ist+1; 
        Prec(ist,1:Ntime)=Rain_Inp(i,j,1:Ntime); ! mm/month
        Temp(ist,1:Ntime)=Temp_Inp(i,j,1:Ntime)  ! C deg
     Enddo
  Enddo
  Deallocate (Xtime)
  print*,"Nst,Ntime",Nst,Ntime
End Subroutine Read_APH_netcdf
!------------------------------------------------------------------------------
Subroutine Write_All_Dr_IDX_to_NetCDF(Fname) !, Nlon, Nlat, Ntime)
   Integer :: ist, jlat, ilon, it, ilen, ix
   !   Integer :: Nlon, Nlat, Ntime
   !Real, Dimension (Nlon, Nlat, Ntime) :: X1, X2, X3, X4, X5, X6
   !Real, Dimension (Ntime) :: Xt
   Real, Allocatable :: Lon(:), Lat(:)
   Character (Len=50) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, t_nm
   Character (Len=*) :: Fname

   print*, Nlon, Nlat, Ntime
   print*, Lon, Lat

   if (.Not.Allocated(X1)) Allocate(X1(Nlon, Nlat, Ntime))
   if (.Not.Allocated(X2)) Allocate(X2(Nlon, Nlat, Ntime))
   if (.Not.Allocated(X3)) Allocate(X3(Nlon, Nlat, Ntime))
   if (.Not.Allocated(X4)) Allocate(X4(Nlon, Nlat, Ntime))
   if (.Not.Allocated(X5)) Allocate(X5(Nlon, Nlat, Ntime))
   if (.Not.Allocated(X6)) Allocate(X6(Nlon, Nlat, Ntime))
   if (.Not.Allocated(Xt)) Allocate(Xt(Ntime))
   Allocate (Lon(Nlon), Lat(Nlat))
   Lon(:)=xlon(:); Lat(:)=xlat(:)

   Do it=1, Ntime;   Xt(it) = it; Enddo
   ist = 0
   ! PDSI(Nst,Ntime):
   Do jlat=1, Nlat
      Do ilon=1, Nlon
         ist = ist +1
         Do it = 1, Ntime
            X1(ilon,jlat,it) = mon_ori_pdsi(ist,it)
            X2(ilon,jlat,it) = mon_ori_wplm(ist,it)
            X3(ilon,jlat,it) = mon_ori_zind(ist,it)
            X4(ilon,jlat,it) = mon_sel_pdsi(ist,it)
            X5(ilon,jlat,it) = mon_sel_wplm(ist,it)
            X6(ilon,jlat,it) = mon_sel_zind(ist,it)
         Enddo
      Enddo
   Enddo

   Write(t_unit,"(A,I4.4,A)") "months since ",YBas1,"-01-01"
   print*, "writing ori_wplm ..."
   Var_nm="ori_pdsi"; Var_ln="Origin PDSI"; var_unit="-"; lon_nm="lon"; lat_nm="lat"; t_nm="month"
   Call wr_ncdf_3d (trim(Fname), Lon, Lat, Xt, Nlon, Nlat, Ntime, X1,trim(Var_nm),trim(t_unit),trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)
   Var_nm="ori_wplm"; Var_ln="Origin wplm"
   Call wr_ncdf_3d_apd (trim(Fname), X2, trim(Var_nm),Nlon, Nlat, Ntime, trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)
   Var_nm="ori_zind"; Var_ln="Origin zind"
   Call wr_ncdf_3d_apd (trim(Fname), X3, trim(Var_nm),Nlon, Nlat, Ntime, trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)
   Var_nm="sel_pdsi"; Var_ln="Self calibrate PDSI"
   Call wr_ncdf_3d_apd (trim(Fname), X4, trim(Var_nm),Nlon, Nlat, Ntime, trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)
   Var_nm="sel_wplm"; Var_ln="Self calibrate wplm"
   Call wr_ncdf_3d_apd (trim(Fname), X5, trim(Var_nm),Nlon, Nlat, Ntime, trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)
   Var_nm="sel_zind"; Var_ln="Self calibrate zind"
   Call wr_ncdf_3d_apd (trim(Fname), X6, trim(Var_nm),Nlon, Nlat, Ntime, trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)

   ! SPI, SPE(Nst,Ntime,NLen)
   Do ilen=1, NLen
      ist=0
      Do jlat=1, Nlat
         Do ilon=1, Nlon
            ist = ist +1
            Do it = 1, Ntime
               X1(ilon,jlat,it) = SPE(ist,it,ilen)
               X2(ilon,jlat,it) = SPI(ist,it,ilen)
            Enddo
         Enddo
      Enddo
      write(Var_nm,"(A,I2.2)") "SPEI_", Len_SPI(ilen)
      write(Var_ln,"(A,I2.2)") "SPEI index in timescale of ", Len_SPI(ilen)

!      if (ilen==1) then
!         Call wr_ncdf_3d (trim(Fname), Lon, Lat, Xt, Nlon, Nlat, Ntime, X1,trim(Var_nm),trim(t_unit),trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)
!      else
         Call wr_ncdf_3d_apd (trim(Fname), X1, trim(Var_nm),Nlon, Nlat, Ntime, trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)
!      endif

      write(Var_nm,"(A,I2.2)") "SPI_", Len_SPI(ilen)
      write(Var_ln,"(A,I2.2)") "SPI index in timescale of ", Len_SPI(ilen)
      Call wr_ncdf_3d_apd (trim(Fname), X2, trim(Var_nm),Nlon, Nlat, Ntime, trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)
   Enddo

   ! Luu Nhiet do va luong mua dau vao:
   Var_ln="Monthly T2m-temperature"
   Call wr_ncdf_3d_apd (trim(Fname), Temp_Inp, "T2m", Nlon, Nlat, Ntime, trim(Var_ln), "Deg.C",trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)
   Var_ln="Monthly Precipitation Count"
   Call wr_ncdf_3d_apd (trim(Fname), Rain_Inp, "Pr", Nlon, Nlat, Ntime, trim(Var_ln), "mm/month",trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)

   ! Luu Mask_Land:
   Var_ln="Land mask: =1: Land, =-99: Sea/water";
   Call wr_ncdf_2d_apd (trim(Fname), Mask, "land", nlon, nlat, trim(Var_ln), "-", trim(lon_nm), trim(lat_nm), amssng)

End Subroutine Write_All_Dr_IDX_to_NetCDF
!------------------------------------------------------------------------------
SUBROUTINE putAWC(dat,lon,lat,nst)
   integer :: nst, ist
   real :: dat(nst), lon(nst), lat(nst)
   integer, parameter :: ncols = 1099, nrows = 549
   integer :: icol, irow, itmp
   !character (len = *), parameter :: fileAWC = "/work/users/thanhnx/Func_Sub/AWC_Global.asc"
   real, parameter :: xcor = -180.0, ycor = -89.9056530001
   real, parameter :: cellsiz = 0.327697, iflag = -9999, oflag = -99
   real :: indat(ncols,nrows), x(ncols), y(nrows)
   ! <-~~~=|==-----------------------------------------------------------<o
   ! <-~~~=|== Read asc AWC file
   !open (1,file = fileAWC,status = "old")
   Open (1, file = trim(AWC_Fname), Status="Old")
      ! <-~~~=| Skip reading 7 lines
      read (1,*) ; read (1,*) ; read (1,*) ; read (1,*)
      read (1,*) ; read (1,*) ; read (1,*)
      ! <-~~~=|
      do irow = 1, nrows
         itmp = nrows - irow + 1
         read (1,*) (indat(icol,itmp),icol = 1, ncols)
      enddo
   close (1)
   do icol = 1, ncols
      x(icol) = xcor + (icol - 1) * cellsiz
   enddo
   do irow = 1, nrows
      y(irow) = ycor + (irow - 1) * cellsiz
   enddo
   ! <-~~~=| Compute for new grid
   do ist = 1, nst
      icol = nrtpos(x,lon(ist),ncols,oflag)
      irow = nrtpos(y,lat(ist),nrows,oflag)
      dat(ist) = indat(icol,irow)*0.0393701
      if (dat(ist) == iflag) then
         dat(ist) = oflag   ! Doi -99.9 --> 99.0 
      endif
   !   print*,ist, dat(ist)
   enddo
return
END SUBROUTINE putAWC
!--------------------------------------------------------------------
integer FUNCTION nrtpos(x,y,n,flag)
   integer :: n, i, itmp
   real :: flag, x(n), y, minx

   minx = 999999
   itmp = 0
   do i = 1, n
      if (x(i) /= flag .and. y /= flag) then
         if (abs(x(i) - y) < minx) then
            minx = abs(x(i) - y)
            itmp = i
         endif
      endif
   enddo

   nrtpos = itmp
   return
END FUNCTION nrtpos
!--------------------------------------------------------------------
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_RW_Station_nc.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.f90" ! Chuong trinh tinh chi so han SPEI

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