!------------------------------------------------------------------------------
Subroutine Write_Time_Series_to_TEXT (X, iunit, VarName)
Integer :: ilen, iyr, iunit
Real :: X(Nyear,NLen)  ! N = Nyear; M = NLen
Character (Len=*) :: VarName
write (iunit,*)
write (iunit,"(A)") "Drought Characteristic: "//trim(VarName)
write (iunit,"(A5,6I8)") "Year",(Len_SPI(ilen), ilen=1,NLen)
Do iyr=1, NYear
   write (iunit,"(I4,1x,6F8.2)") iyr+YBas1-1, (X(iyr,ilen),ilen=1,Nlen)
Enddo
End Subroutine Write_Time_Series_to_TEXT
!------------------------------------------------------------------------------
Subroutine Write_2D_to_NC (X, Grads_File, Field, Append)
! Chuyen mot vecto chua cac diem luoi ve mang 2 chieu va ghi ra file Grads
Integer :: i,j,i1, ii
Real :: X(Nst_Max), Lon0, Lat0, DelX !, Xlon(Nlon), Xlat(Nlat)
Real :: Y(Nlon, Nlat), flag
Character (Len=*) :: Grads_File, Field  ! Grads file, Characteristic, Drought scale
Logical :: Append
flag = 1.E+38 !-99999.
i1 = 0
Do j=1, Nlat
   Do i=1, Nlon
      i1 = i1 + 1
      Y(i,j) = X(i1)
   Enddo
Enddo
Where (Y==-99) Y=flag
!Xlon(1:Nlon) = Lon(1:Nlon); Xlat(1:Nlat) =  Lat(1:Nlat)
if (Append) then
   Call wr_ncdf_2d_apd (trim(Grads_File), Y, trim(Field), Nlon,Nlat,flag)
else
   call wr_ncdf_2d (trim(Grads_File),Xlon,Xlat,Nlon,Nlat,Y,trim(Field),flag)
endif
End Subroutine Write_2D_to_NC
!------------------------------------------------------------------------------
Subroutine Write_SPI_To_TEXT (Fname)
Character (Len=*) :: Fname
Integer :: ifirst, last, i, iy, im, ist
Open (3, File=trim(Fname), Status="Unknown")
Do ist=1, Nst
!   write (3,"(A)") Station(ist)  
!   write (3, "(A4,A3,4I7)") "Yr","Mn", (Len_SPI(i),i=1,NLen)
!  Skip leading missings
   ifirst = 0
   Do
      ifirst = ifirst + 1
      if (spi(ist,ifirst,1) /= amssng) exit
   Enddo
!  Skip trailing missngs
   last = maxyrs*Nmon_week + 1
   Do
      last = last - 1
      if (spi(ist,last,1) /= amssng) exit
   Enddo
Enddo
!  Ouput SPI's
Do j=1, NLen
   write (3, "(A,I2.2,A)") "      Drought time scale: ", Len_SPI(j), " mon"
   write (3, "(A8,100(1x,A7))") "Mon\Stn",(trim(Station(ist)), ist=1, Nst), (trim(Station(ist)), ist=1, Nst)
   Do i = ifirst, last
      iy = (i-1)/Nmon_week + YBas1 ! ibegyr
      im = mod (i-1,Nmon_week)+1
      write(3, "(i4,i2.2,2x,100f8.2)") iy, im, (spi(ist,i,j), ist=1, Nst), (DatBL(i,ist),ist=1, Nst) ! j=1,NLen)
   Enddo
Enddo
End Subroutine
!------------------------------------------------------------------------------
Subroutine Write_Frequency_TEXT (Res_Dir, FScale)
Character (Len=255) :: FTime, FSpace
Character (Len=*) :: Res_Dir, FScale
Integer :: ilen, ifirst, last, i, iy, im, ist

! SPI
Do ilen=1, NLen
   write (FTime,"(A,I2.2,A)") trim(Res_Dir)//"SPI_Time_Freq"//trim(FScale),Len_SPI(ilen),".txt"
   write (FSpace,"(A,I2.2,A)") trim(Res_Dir)//"SPI_Spac_Freq"//trim(FScale),Len_SPI(ilen),".txt"
   Open (3, File=trim(FTime), Status="Unknown")
   Open (4, File=trim(FSpace), Status="Unknown")
   write (4,"(A21,11A9)") "Station     ",(SPI_Lev(i),i=1,NThre_SPI+1)
   write (3,"(A4,A3,11A9)") "Yr","Mn",(SPI_Lev(i),i=1,NThre_SPI+1)
   Do ist=1, Nst
      write (4,"(A21,11F9.3)") Station(ist), (SPI_Spac_Freq(ist,i,ilen),i=1,NThre_SPI+1)
   Enddo

   Do it = 1, NTime
      iy = (it-1)/Nmon_week + YBas1; im = mod (it-1,Nmon_week)+1
      write(3,"(i4,i3,11f9.3)") iy, im, (SPI_Time_Freq(it,i,ilen),i=1,NThre_SPI+1)
   Enddo
   Close (3); Close(4)
Enddo
End Subroutine
!------------------------------------------------------------------------------
Subroutine Plot_Hovm_Freq (FScale)
Real, Allocatable, Dimension (:) :: X
Integer :: i, ist, ithre, it, Nx, Ny, ilen
Character (Len=100) :: GrSt
Character (Len=*) :: FScale

!Time SPI  
Nx=Ntime; Ny=NThre_SPI+1; Allocate (X(Nx*Ny))
Do ilen=1, NLen
   write (GrSt,"(A,I2.2)") "SPI_T_Freq"//trim(FScale),Len_SPI(ilen); i=0
   Do ithre=1, Ny;   Do it=1, Nx;      i=i+1; X(i)=SPI_Time_Freq(it,ithre,ilen); Enddo; Enddo
   Call Hovmoller_Freq (X, Nx, Ny, trim(GrSt), "./4Grads/", trim(GrSt))
Enddo
Deallocate (X)
! Spacial SPI
Nx=Nst; Ny=NThre_SPI+1; Allocate (X(Nx*Ny))
Do ilen=1, NLen
   write (GrSt,"(A,I2.2)") "SPI_S_Freq"//trim(FScale),Len_SPI(ilen); i=0
   Do ithre=1, Ny;   Do ist=1, Nx;      i=i+1; X(i)=SPI_Spac_Freq(ist,ithre,ilen); Enddo; Enddo
   Call Hovmoller_Freq (X, Nx, Ny, trim(GrSt), "./4Grads/", trim(GrSt))
Enddo
Deallocate (X)
End Subroutine
!------------------------------------------------------------------------------
Subroutine Plot_Hovm_SPI
Integer :: it1, it2, i
Character (Len=50) :: GrFile
! Plot SPI for all Len_SPI (1,3,6,12)
  it1=1; it2=120   ! 10 nam dau
  Do i=1, NLen
     write (GrFile, "(A,I2.2,A)") "SPI"//trim(FScale),Len_SPI(i),"_1"
     Call Hovmoller_BL (SPI(:,:,i), it1, it2, trim(GrFile),"./4Grads/",trim(GrFile))
  Enddo
  it1=121; it2=240   ! 10 nam tiep
  Do i=1, NLen
     write (GrFile, "(A,I2.2,A)") "SPI"//trim(FScale),Len_SPI(i),"_2"
     Call Hovmoller_BL (SPI(:,:,i), it1, it2, trim(GrFile),"./4Grads/",trim(GrFile))
  Enddo
  it1=241; it2=372   ! 10 nam tiep
  Do i=1, NLen
     write (GrFile, "(A,I2.2,A)") "SPI"//trim(FScale),Len_SPI(i),"_3"
     Call Hovmoller_BL (SPI(:,:,i), it1, it2, trim(GrFile),"./4Grads/",trim(GrFile))
  Enddo
End Subroutine
!------------------------------------------------------------------------------
Subroutine Cal_Drought_Freq
Real, Allocatable, Dimension (:) :: X, Bin, Freq
Integer :: ist, it, ilen

SPI_Time_Freq(:,:,:) = amssng; SPI_Spac_Freq(:,:,:) = amssng
! Time series frequency
Allocate (X(1:Nst), Bin(1:NThre_SPI), Freq(1:NThre_SPI+1))
Do ilen=1, Nlen
   Do it=1, Ntime
      X(1:Nst) = SPI(1:Nst,it,ilen); Bin(1:NThre_SPI)=SPI_Thre(1:NThre_SPI)
      Call Frequency (X, Nst, Bin, Freq, NThre_SPI)
      SPI_Time_Freq(it,1:NThre_SPI+1,ilen) = Freq(1:NThre_SPI+1)
   Enddo
Enddo
Deallocate (X,Bin,Freq)
! Spacial frequency
Allocate (X(1:Ntime), Bin(1:NThre_SPI), Freq(1:NThre_SPI+1))
Do ilen=1, Nlen
   Do ist=1, Nst
      X(1:Ntime) = SPI(ist,1:Ntime,ilen); Bin(1:NThre_SPI)=SPI_Thre(1:NThre_SPI)
      Call Frequency (X, Ntime, Bin, Freq, NThre_SPI)
      SPI_Spac_Freq(ist,1:NThre_SPI+1,ilen) = Freq(1:NThre_SPI+1)
   Enddo
Enddo
Deallocate (X,Bin,Freq)

End Subroutine
!------------------------------------------------------------------------------
Subroutine Frequency (X, N, Bin, Freq, M)
Integer :: N, M, i, j, ic
Real :: X(N), Bin(M), Freq(M+1), tmp
ic=0
Do i=1, N
   if (X(i) /= amssng) ic=ic+1
Enddo
tmp=ic
Do j=1, M
   if (j==1) then
      ic=0
      Do i=1, N
         if (X(i) /= amssng .AND. X(i) <= Bin(j)) ic=ic+1
      Enddo
   else
      ic=0
      Do i=1, N
         if (X(i) > Bin(j-1) .AND. X(i) <= Bin(j)) ic=ic+1
      Enddo
   endif
   if (tmp==0) then; Freq(j) = amssng; else; Freq(j) = real(ic)/tmp; endif
Enddo
ic=0
Do i=1, N
   if (X(i) > Bin(M)) ic=ic+1
enddo
if (tmp==0) then; Freq(M+1) = amssng; else; Freq(M+1) = real(ic)/tmp; endif
End Subroutine
!------------------------------------------------------------------------------
Subroutine Hovmoller_BL (X, it1, it2, Grads_File, Out_Dir, Field)
Integer :: it, it1, it2, ist, Nt, irec
Character (Len=*) :: Grads_File, Out_Dir, Field
Real :: X(Nst_Max,NYear*NWeek), Y(Ntime,Nst)
 Open (33,file=trim(Out_Dir)//trim(Grads_File)//trim(FScale)//'.dat',FORM='UNFORMATTED', &
          access="direct",recl=4,status="unknown") !*Nt*Nst,status="unknown")
! Where (X > 50) X=-99.; Where (X /= -99 .and. X <-50) X=-99.
 Nt = it2-it1+1; irec = 0
 Do ist=1, Nst
!    Y(1:Ntime,ist) = X(ist,1:Ntime)
    Y(1:Nt,ist) = X(ist,it1:it2)
 Enddo
 Do it=1, Nt
    Do ist=1,Nst
       irec=irec+1; write (33,rec=irec) Y(it,ist)
!       write (33,rec=1) Y
    Enddo
 Enddo
 Close (33)
 Open (4, file=trim(Out_Dir)//trim(Grads_File)//trim(FScale)//".ctl")
 write (4,"(A)") "dset ^"//trim(Grads_File)//trim(FScale)//".dat"
 write (4,"(A)") "title Spatial-Temporal Distribution of "//trim(Field)
 write (4,"(A,F10.3)") "undef -99.0"
 write (4,"(A,I4,A,I4,A)") "xdef ",Nst," linear ",1," 1"
 write (4,"(A,I4,A)") "ydef ",Nt," linear 1 1"
! write (4,"(A,I4,A)") "ydef ",1," linear 1 1"
! write (4,"(A)") "options yrev"
 write (4,"(A)") "zdef 1 linear 1000 1"
! write (4,"(A,I4,A,I4,A)") "tdef ",Ntime," linear 00z01Jan",YBas1," 1mon"
 write (4,"(A,I4,A,I4,A)") "tdef ",1," linear 00z01Jan",YBas1," 1mon"
 write (4,"(A,I3)") "vars 1"
 write (4,"(A,A,A)") trim(Field)," 0 99 "," Variable "
 write (4,"(A)") "endvars"
 Close(4)
End Subroutine Hovmoller_BL
!------------------------------------------------------------------------------
Subroutine Hovmoller_Freq (X, Nx, Ny, Grads_File, Out_Dir, Field)
Integer :: Nx, Ny, i, irec
Character (Len=*) :: Grads_File, Out_Dir, Field
Real :: X(Nx*Ny)
 Open (33,file=trim(Out_Dir)//trim(Grads_File)//trim(FScale)//'Freq.dat',FORM='UNFORMATTED', &
          access="direct",recl=4,status="unknown") !*Nt*Nst,status="unknown")
! Where (X > 50) X=-99.; Where (X /= -99 .and. X <-50) X=-99.
 irec = 0
 Do i=1, Nx*Ny
    irec=irec+1; write (33,rec=irec) X(i)
 Enddo
 Close (33)
 Open (4, file=trim(Out_Dir)//trim(Grads_File)//trim(FScale)//"Freq.ctl")
 write (4,"(A)") "dset ^"//trim(Grads_File)//trim(FScale)//"Freq.dat"
 write (4,"(A)") "title Frequency Distribution of "//trim(Field)
 write (4,"(A,F10.3)") "undef -99.0"
 write (4,"(A,I4,A,I4,A)") "xdef ",Nx," linear ",1," 1"
 write (4,"(A,I4,A)") "ydef ",Ny," linear 1 1"
! write (4,"(A,I4,A)") "ydef ",1," linear 1 1"
! write (4,"(A)") "options yrev"
 write (4,"(A)") "zdef 1 linear 1000 1"
! write (4,"(A,I4,A,I4,A)") "tdef ",Ntime," linear 00z01Jan",YBas1," 1mon"
 write (4,"(A,I4,A,I4,A)") "tdef ",1," linear 00z01Jan",YBas1," 1mon"
 write (4,"(A,I3)") "vars 1"
 write (4,"(A,A,A)") trim(Field)," 0 99 "," Variable "
 write (4,"(A)") "endvars"
 Close(4)
End Subroutine Hovmoller_Freq
!------------------------------------------------------------------------------