! 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 = 1986, YBas2 = 2005, 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 

Character (Len=255) :: Fname, 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 (:,:,:) :: Prec
!------------------------------------------------------------------------------
Integer, Parameter :: NFile = 30  ! 6 models + 6 models * 2 Scenarios * 2 Future Periods
Character (Len=100) :: Fname_inp(NFile), Fname_Out(NFile)
Namelist /Filenames/ Fname_inp, Fname_Out
Integer :: ifile
!------------------------------------------------------------------------------

!DatBL(:,:) = amssng ! (Ntime, Nst_Max)
!SPI(:,:,:) = amssng
!------------------------------------------------------------------------------
Res_Dir = "../RES/SPI_Res/"

!Call Read_data_netcdf

! Lay thong tin cac vung khi hau: B1=2, B2=3,..., N3=8
Call rd_ncdf_2d ("/home/tanpv/SEACLID/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

Open (1, file = "File_List.txt"); Read (1, Filenames); Close(1)
Do ifile=1, NFile
   DatBL(:,:) = amssng ! (Ntime, Nst_Max)
   SPI(:,:,:) = amssng
   Fname = Fname_inp(ifile)
   Call Read_data_netcdf(Fname,ifile)
   ! 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 "(12F12.3)", (SPI(ist,it,1),it=1,Ntime)
   !   pause
   !Enddo
   ! Print SPI to file
   Fname=trim(Fname_Out(ifile))
   Call Write_SPI_to_NetCDF (Fname)
   print*,"Done: Writ to ",Fname
Enddo
print*,"Program successfulll terminated !!!"

!==============================================================================
CONTAINS

!------------------------------------------------------------------------------
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

t_unit="months since 1986-01-15"
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), Prec, "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(FileName,ifile)
Character (Len=*) :: FileName 
Character (Len=255) :: VarName
Real, Allocatable :: Xlon(:), Xlat(:), XTime(:) !, X(:,:,:)
Real, allocatable, dimension (:,:,:) :: Precip
!Real, Allocatable :: XTime(:), X(:,:,:)
Integer :: Nlon1, Nlat1, Nt, i, j, ist, it, iyr, imon, ifile

!  FileName="../Data/APH_hist_allyear_40y.nc"
  VarName ="pr"
  Call rd_ncdf_3d (trim(FileName), Xlon, Xlat, XTime, Nlon1, Nlat1, Nt, Precip, trim(VarName))
  if (.not.Allocated(XXlon)) Allocate (XXlon(Nlon1)); if (.not.Allocated(XXlat)) Allocate (XXlat(Nlat1))
  if (.not.Allocated(Prec)) Allocate (Prec(Nlon1,Nlat1,Nt))
  ! Allocate (XXlon(Nlon1), XXlat(Nlat1))
  XXLon(:) = Xlon(:); XXLat(:) = Xlat(:); Prec(:,:,:) = Precip(:,:,:)

  print*,"In Read_data_netcdf: Nlon, Nlat, Nt=", Nlon1, Nlat1, Nt
  print*,Xlon, Xlat
  Where (Precip < 0) Precip = amssng; 
  ist = 0
  Do j=1, Nlat1
     Do i=1, Nlon1
        ist=ist+1;       !write(Station(ist),"(A,I4.4)") "Point-",ist
        Lon(ist)=Xlon(i); Lat(ist)=Xlat(j);
     !   VN_Land(ist)=vietnam(i,j)
        if (ifile<=6) then; iyr=1985; elseif (ifile<=18) then; iyr=2045; elseif (ifile<=30) then; iyr=2079; endif
     !   iyr=1985
        Do it=1, Nt
           imon=mod(it,NMonth); if (imon==1) iyr=iyr+1; if (imon==0) imon=12
           if (Precip(i,j,it)/=amssng) DatBL(it,ist) = Precip(i,j,it)*DayOfMon(imon,iyr)
!           print*,iyr,imon,DayOfMon(imon,iyr)
        Enddo
!        pause
!        print "(12F6.1)",DatBL(:,ist)
     Enddo
  Enddo
  Nst = ist; Nlon = Nlon1; Nlat = Nlat1
  print*,"Nst,Ntime",Nst,Nt
!  pause
End Subroutine Read_data_netcdf
!------------------------------------------------------------------------------
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 "Sub_SPI.f90" ! Chuong trinh tinh chi so han SPI
Include "sub_wr_ncdf.F90"
Include "sub_rd_ncdf.F90"
!Include "../Libs/sub_wr_ncdf_Tan.F90"
!==============================================================================
END   ! Program