! Calculate SPI and write to NetCDF file for further use
! Using Monthly data only
use netcdf
Implicit None
Integer, Parameter :: Nst_Max = 1000   ! So Grid Cells hoac so tram
Integer, Parameter :: NMonth = 12, NLen = 4
Integer, Parameter :: Len_SPI(NLen) = (/ 1, 3, 6, 12 /) ! Cac dot han 1,2,3,4,12 thang/tuan theo SPI 
Integer, Parameter :: YBas1 = 1961, YBas2 = 2014, NYear = YBas2-YBas1+1, maxyrs = NYear

! For testing with origin data:
!Integer, Parameter :: YBas1 = 1950, YBas2 = 1992, NYear = YBas2-YBas1+1, maxyrs = NYear

Integer, Parameter :: Ntime = NYear*NMonth

! So lieu Khi tuong dau vao
Real, Dimension (Ntime) :: Rain     ! Meteo Input Data for Baseline
! SPI
Real, Parameter :: amssng = -99.0
Real, Dimension (NMonth) :: beta_spi, gamm, pzero
Real, Dimension (Nst_Max, Ntime, NLen) :: SPI
!Thong tin tram
Integer :: NSt, idx, Nlon, Nlat
Real, Dimension (Nst_Max) :: Lat, Lon 
Real, Allocatable, Dimension (:) :: XXLat, XXLon
Real, Dimension (Ntime, Nst_Max) :: DatBL, DatTMP 

Character (Len=255) :: Fname, BL_Dir, Res_Dir, FScale, VarName
Integer :: ist, it, j
Logical :: status

Integer :: nxlon, nxlat, ilon, jlat
Real, allocatable, dimension (:)   :: loninfo, latinfo, Xlon, Xlat
Real, allocatable, dimension (:,:) :: vietnam ! thong tin cac vung khi hau: B1=2, B2=3,..., N3=8
Real, allocatable, dimension (:,:,:) :: Precip
Real, Allocatable :: XTime(:), Preci(:,:), XH(:), Xreg(:)
Character (Len=20), allocatable, dimension (:) :: Stname
!------------------------------------------------------------------------------
DatBL(:,:) = amssng ! (Ntime, Nst_Max)
DatTMP(:,:) = amssng
SPI(:,:,:) = amssng
!------------------------------------------------------------------------------
BL_Dir  = "/home/tanpv/SEACLID_New/Data/"
Res_Dir = "SPI_RES/"

!Call Read_data_Baseline

Call Read_data_netcdf

!Call Read_Data_test_SPI   ! Test with origin data

! Lay thong tin cac vung khi hau: B1=2, B2=3,..., N3=8
!Call rd_ncdf_2d ("/home/tanpv/SEACLID_New/MASK/VN.grid_0.25.nc",loninfo,latinfo,nxlon,nxlat,vietnam,"all")
!Where (vietnam==2) vietnam=1; Where (vietnam==3) vietnam=2; Where (vietnam==4) vietnam=3;
!Where (vietnam==5) vietnam=4; Where (vietnam==6) vietnam=5; Where (vietnam==7) vietnam=6;
!Where (vietnam==8) vietnam=7   ! Gan lai cac vung KH: B1=1, B2=2, B3=3, B4=4, N1=5, N2=6, N3=7
  
! SPI
Print*,"Calculating, Writing and Ploting SPI ..."
Call Cal_SPI
Do ist=1, Nst
!   if (Empty_Ar1 (SPI(ist,1:Ntime,1), Ntime)) cycle
!   print"(12F7.3)", (SPI(ist,it,1),it=1,Ntime)
!   pause
Enddo
! Print SPI to file

Fname=trim(Res_Dir)//"SPI_Stn_B1-3"
Call Write_SPI_to_NetCDF (Fname)

SPI(:,:,:) = amssng
Call Cal_SPEI

Fname=trim(Res_Dir)//"SPEI_Stn_B1-3"
Call Write_SPI_to_NetCDF (Fname)

! Call Write_Results_test
print*,"Done ..."
!==============================================================================
CONTAINS
! Include "/home/tanpv/SEACLID_New/SPEI/OLD/Read_Data_test.f90" ! For test only !!
!------------------------------------------------------------------------------
Subroutine Write_SPI_to_NetCDF (Fname)
! Real, Dimension (Nst_Max, Ntime, NLen) :: SPI
! wr_ncdf_3d_apd (filenm,var,var_nm,nlon,nlat,nt,flag)
! wr_ncdf_3d (filenm,lon,lat,t,nlon,nlat,nt,var,var_nm,flag)
Integer :: jlat, ilon, it, ilen
Real, Allocatable, Dimension (:,:) :: X
Real, Allocatable, Dimension (:)   :: XStn
!Real :: Xlon(Nlon), Xlat(Nlat)
Real, Dimension (Ntime) :: t
Logical :: Append
Character (Len=*) :: Fname
Character (Len=30) :: t_unit

t_unit="months since 1961-01-01"
print*,"Xlon=",Xlon
Allocate (X(Nst,Ntime),XStn(Nst))
Do it=1, Ntime;   t(it) = it-1 ; Enddo;   Do ist=1, Nst; XStn(ist) = ist; Enddo

Append = .False.
Do ilen=1, Nlen
   Do ist=1, Nst
      Do it=1, Ntime
         X(ist,it) = SPI(ist,it,ilen)
     !       print*, X(ilon,jlat,it)
      Enddo
   Enddo
   write(VarName,"(A,I2.2)") "SPI_", Len_SPI(ilen)
   if (Append) then
      Call wr_ncdf_2d_apd (trim(Fname)//".nc", X, trim(VarName), Nst, Ntime, "SPI", "-", "Station", "Time", amssng)
      Append = .True.
   else
      Call wr_ncdf_2d (trim(Fname)//".nc", XStn, t, Nst, Ntime, X, trim(VarName), "SPI", "-", "Station", "Time", amssng)
      Append = .True.
   endif
Enddo
print*,"writing SPI done"

Call wr_ncdf_1d_apd (trim(Fname)//".nc", Xlon, "lon", Nst, "longitude", "deg. E", "Station", amssng)
print*,"lon>>>"
Call wr_ncdf_1d_apd (trim(Fname)//".nc", Xlat, "lat", Nst, "latitude", "deg. N", "Station", amssng)
Call wr_ncdf_1d_apd (trim(Fname)//".nc", Xreg, "reg", Nst, "region", "-", "Station", amssng)
Call wr_ncdf_1d_apd_char (trim(Fname)//".nc",Stname,"stname", Nst, "Station name","-", "Station", amssng)
Deallocate (X, XStn)

End Subroutine Write_SPI_to_NetCDF
!------------------------------------------------------------------------------
Subroutine Read_data_netcdf
Character (Len=255) :: FileName, VarName
!Real, Allocatable :: XTime(:), Preci(:,:), XH(:), Xreg(:)
!Xlon(:), Xlat(:), XTime(:), Preci(:,:), XH(:), Xreg(:)
Real, Allocatable :: Stnm(:),years(:),mons(:), X(:,:,:), Y(:,:,:)
Integer :: Nlon1, Nlat1, Nt, i, j, ist, it, Nstn, Nyrs, Nmn
Integer :: iyr, imon

  !FileName="/home/tanpv/SEACLID_New/Data/OBS_Stn/R_T2m_1980-2014_mon_Inland_No_Rain_Center.nc"  ! Cac tram tren dat lien
  ! Tinh cho khu vuc Bac Bo (B1-B3)
  FileName="/home/tanpv/OBS_Stn/Monthly/NC/R_T2m_1961-2019_mon_Tram_KT_B1-B3_3D.nc"

!   " R_T2m_1980-2014_mon.nc"
  VarName ="R"
  Call rd_ncdf_3d (FileName,Stnm,years,mons,Nstn,Nyrs,Nmn,X,trim(VarName))
  Deallocate(Stnm,years,mons)
  VarName ="T2m"
  Call rd_ncdf_3d (FileName,Stnm,years,mons,Nstn,Nyrs,Nmn,Y,trim(VarName))
  Deallocate(Stnm)
  Call rd_ncdf_1d (FileName,Stnm, Nstn, Xlon, "Lon")
  Deallocate(Stnm)
  Call rd_ncdf_1d (FileName,Stnm, Nstn, Xlat, "Lat")
  Deallocate(Stnm)
  Call rd_ncdf_1d (FileName,Stnm, Nstn, XH, "H")
  Deallocate(Stnm)
  Call rd_ncdf_1d (FileName,Stnm, Nstn, Xreg, "Reg")
  Deallocate(Stnm)
  Call rd_ncdf_1d_Char (FileName, Stnm, Nstn, Stname, "Stname")
!  Call rd_ncdf_3d (trim(FileName),Xlon,Xlat,XTime, Nlon1, Nlat1, Nt, Precip, trim(VarName))
!  Call read_netcdf_station(trim(FileName),Nst, Nt ,Preci,VarName,XTime, Xlon,Xlat,Xh,Xreg)

!  Allocate (XXlon(Nlon1), XXlat(Nlat1))
!  XXLon(:) = Xlon(:); XXLat(:) = Xlat(:)

!  print*,"In Read_data_netcdf: NStn, Nt=", Nst, Nt
!  print*,Xlon; print*, Xlat; print*,Xh; print*,Xreg
!  Where (Precip < 0) Precip = amssng;
  print*,"Nyrs,Nmn=",Nyrs,Nmn;pause
  Nst=Nstn 
  Do ist=1, Nst
     !DatBL(1:Ntime,ist) = Preci(ist,1:Ntime)
     Lon(ist)=Xlon(ist); Lat(ist)=Xlat(ist);
     it=0
     Do iyr=1,Nyrs
        Do imon=1,Nmn
           it=it+1
           DatBL(it,ist)=X(ist,iyr,imon)
           DatTMP(it,ist)=Y(ist,iyr,imon)
        Enddo
     Enddo
!     print "(12F7.1)", DatBL(1:Ntime,ist); pause
  Enddo
  Nt=it
  print*,"Nst,Nt,Ntime",Nst,Nt,Ntime
  pause
End Subroutine Read_data_netcdf
!------------------------------------------------------------------------------
Include "/home/tanpv/Mylibs/Sub_SPI.f90" ! Chuong trinh tinh chi so han SPI
Include "/home/tanpv/Mylibs/Sub_SPEI.f90"
Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
Include "/home/tanpv/Mylibs/Sub_RW_Station_nc.f90"
!==============================================================================
END   ! Program