! Calculate SPI and write to NetCDF file for further use
! Using Monthly data only
use netcdf
Implicit None
Integer, Parameter :: Nst_Max = 5000   ! 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 = 1980, YBas2 = 2014, NYear = YBas2-YBas1+1, maxyrs = NYear
Integer, Parameter :: Ntime = NYear*NMonth

! Integer, Parameter :: Ntime = 1296  ! For testing with origin data

! So lieu Khi tuong dau vao
Real, Dimension (0:Ntime-1) :: Rain, T2m     ! 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(:,:), Temp(:,:), XH(:), Xreg(:)
!------------------------------------------------------------------------------
DatBL(:,:) = amssng; DatTMP(:,:) = amssng ! (Ntime, Nst_Max)
SPI(:,:,:) = amssng
!------------------------------------------------------------------------------
BL_Dir  = "/home/tanpv/SEACLID_New/Data/"
Res_Dir = "../RES/SPI_Res/"

!Call Read_data_Baseline
 Fname = "/home/tanpv/OBS_Stn/Monthly/NC/R_T2m_1980-2014_mon_Inland.nc"
 Call Read_data_netcdf (Fname)

!T Call Read_Data_test   ! Test with origin data

!Do it=1, Ntime
!   write (*,"(I4,50F8.1)") it, (DatTMP(it,ist),ist=1,10)
!   if (mod(it,50)==0) pause
!Enddo

! 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_SPEI
!Do ist=1, Nst
!   if (Empty_Ar1 (SPI(ist,1:Ntime,1), Ntime)) cycle
!   print"(12F7.1)", (SPI(ist,it,1),it=1,Ntime)
!   pause
!Enddo

! Print SPI to file
 Fname=trim(Res_Dir)//"SPEI_Stn_1980-2014_mon_Inland.nc"
 Call Write_SPI_to_NetCDF (Fname)

!T Call Write_Results_test

print*,"Done ..."
!==============================================================================
CONTAINS
Include "OLD/Read_Data_test.f90"
!------------------------------------------------------------------------------
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, ix
!Real, Dimension (Nlon, Nlat, Ntime) :: X
!Real :: Xlon(Nlon), Xlat(Nlat)
!Real, Dimension (Ntime) :: t
!Logical :: Append
!Character (Len=*) :: Fname
!Character (Len=30) :: t_unit

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 1980-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), X, trim(VarName), Nst, Ntime, "SPI", "-", "Station", "Time", amssng)
      Append = .True.
   else
      Call wr_ncdf_2d (trim(Fname), 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), Xlon, "lon", Nst, "longitude", "deg. E", "Station", amssng)
print*,"lon>>>"
Call wr_ncdf_1d_apd (trim(Fname), Xlat, "lat", Nst, "latitude", "deg. N", "Station", amssng)
Call wr_ncdf_1d_apd (trim(Fname), Xreg, "reg", Nst, "region", "-", "Station", amssng)
Deallocate (X, XStn)


!Do it=1, Ntime
!   t(it) = it-1
!Enddo
!Xlon(1:Nlon) = Lon(1:Nlon); Xlat(1:Nlat) = Lat(1:Nlat)
!Append = .False.
!Do ilen=1, Nlen
!   ix=0
!   Do jlat=1, Nlat
!      Do ilon=1, Nlon
!         ix = ix + 1
!         Do it=1, Ntime
!            X(ilon,jlat,it) = SPI(ix,it,ilen)
!     !       print*, X(ilon,jlat,it)
!         Enddo
!      Enddo
!   Enddo
!   write(VarName,"(A,I2.2)") "SPI_", Len_SPI(ilen)
!   if (Append) then
!      Call wr_ncdf_3d_apd (trim(Fname), X, trim(VarName), Nlon, Nlat, Ntime, amssng)
!      Append = .True.
!   else
!      Call wr_ncdf_3d (trim(Fname), XXLon, XXLat, t, Nlon, Nlat, Ntime, X, trim(VarName), trim(t_unit), amssng)
!      Append = .True.
!   endif
!Enddo
!Call wr_ncdf_3d_apd (trim(Fname), Precip, "Precip", Nlon, Nlat, Ntime, amssng)
!Call wr_ncdf_2d_apd (trim(Fname), vietnam, "Subreg", nlon, nlat, amssng)
End Subroutine Write_SPI_to_NetCDF
!------------------------------------------------------------------------------
Subroutine Read_data_netcdf (Fname)
Character (Len=*) :: Fname
Character (Len=255) :: VarName
!Real, Allocatable :: Xlon(:), Xlat(:), XTime(:) !, X(:,:,:)
!Real, Allocatable :: XTime(:), X(:,:,:)
Integer :: Nlon1, Nlat1, Nt, i, j, ist, it

!  FileName="/home/tanpv/SEACLID_New/Data/OBS_Stn/R_T2m_1980-2014_mon_Inland_No_Rain_Center.nc"
!" R_T2m_1980-2014_mon.nc"
  VarName ="R"
  Call read_netcdf_station(trim(FName),Nst, Nt ,Preci, VarName, XTime, Xlon,Xlat,Xh,Xreg)

  Deallocate (Xlon,Xlat,Xh,Xreg)
  VarName ="T2m"
  Call read_netcdf_station(trim(FName),Nst, Nt ,Temp, VarName, XTime, Xlon,Xlat,Xh,Xreg)

  Where (IsNaN(Preci)); Preci = amssng; Endwhere; Where (IsNaN(Temp)); Temp = amssng; Endwhere 
  Call FillUp_Miss_Data (Temp, Ntime, Nst); Call FillUp_Miss_Data (Preci, Ntime, Nst)

  Do ist=1, Nst
     DatBL(1:Ntime,ist) = Preci(ist,1:Ntime); DatTMP(1:Ntime,ist) = Temp(ist,1:Ntime)
     Lon(ist)=Xlon(ist); Lat(ist)=Xlat(ist);
!     print "(12F7.1)", DatTMP(1:Ntime,ist); pause
  Enddo
  print*,"Nst,Ntime",Nst,Nt
End Subroutine Read_data_netcdf
!------------------------------------------------------------------------------
Subroutine FillUp_Miss_Data (X, Ntime, Nst)
Integer :: ist, it, imon, it1, iyr, Nst, Ntime
Real :: X(Nst,Ntime), Mean(NMonth)
Do ist=1, Nst
   Do imon=1, Nmonth
      Mean(imon)=0;iyr=0
      Do it=imon, Ntime, Nmonth
         if (X(ist,it)/=amssng) then; Mean(imon)=Mean(imon)+X(ist,it); iyr=iyr+1;endif
      Enddo
      Mean(imon)=Mean(imon)/real(iyr)
      Do it=imon, Ntime, Nmonth
         if (X(ist,it)==amssng) then; print*,X(ist,it); X(ist,it)=Mean(imon); print*,X(ist,it); pause;endif
      Enddo
   Enddo
Enddo
End Subroutine FillUp_Miss_Data
!------------------------------------------------------------------------------
Include "Sub_SPEI.f90" ! Chuong trinh tinh chi so han SPI
Include "../Libs/sub_wr_ncdf_Tan.F90"
Include "../Libs/sub_rd_ncdf.F90"
Include "../Libs/Sub_RW_Station_nc.f90" 
!==============================================================================
END   ! Program