! Chuong trinh tinh chi so han PDSI
! Hai buoc tinh chinh:
! 1. Tinh cac he so cho thoi ky Baseline (thuong chon 20 nam)
! 2. Tinh cac chi so Z va Vpml de du bao cho tuong lai
! Qua trinh tinh thuc hien cho tung tram. Co the tinh theo thang hoac tuan
! So lieu dau vao bat buoc phai co (Trung binh/tong luong Thang hoac Tuan 10 ngay):
! H, Lat, Lon, AWC, Rain, Tx, Tm
! Day la phien ban da chinh sua theo cach moi 
Use NetCDF
Implicit None
Character (Len=5)  :: Decade
Integer, Parameter :: Nst_Max = 5000   ! So Grid Cells hoac so tram
Integer, Parameter :: NMonth = 12, NWeek = 52, NLen = 4
Real,    Parameter :: Nguong = -1.5    ! Nguong han nang 
Integer, Parameter :: NThre_SPI = 10
Real, Parameter :: SPI_Thre(NThre_SPI) = (/-2.,-1.6,-1.3,-0.8,-0.5,0.5,0.8,1.3,1.6, 2./)
Character (Len=7), Parameter :: SPI_Lev(NThre_SPI+1) = (/"Exc_Dry","Ext_Dry","Sev_Dry","Mod_Dry","Abn_Dry", &
                                            "Nea_Nor","Abn_Moi","Mod_Moi","Ver_Moi","Ext_Moi","Exc_Moi"/)
Integer, Parameter :: NdayInMonth (NMonth) = (/31,28,31,30,31,30,31,31,30,31,30,31/)
!SPI
Integer, Parameter :: Len_SPI(NLen) = (/ 1, 3, 6, 12 /) ! Cac dot han 1,2,3,4,12 thang/tuan theo SPI 
Integer, Parameter :: YBas1 = 1956, YBas2 = 2005, NYear = YBas2-YBas1+1, maxyrs = NYear
!Integer, Parameter :: YBas1 = 1950, YBas2 = 1992, NYear = YBas2-YBas1+1, maxyrs = NYear 
                       ! SL 20 nam chuan de tao SL tron ==> YBas2-YBas1+1=20
Character (Len=21) :: Station(Nst_Max), Stn_Fcs(Nst_Max) ! Danh sach tram ban dau va danh sach tram Du bao
! Drought Frequency
Real, Dimension (NYear*NWeek,NThre_SPI+1,NLen) :: SPI_Time_Freq
Real, Dimension (Nst_Max,NThre_SPI+1,NLen) :: SPI_Spac_Freq
!
! So lieu Khi tuong dau vao
Real, Dimension (NYear*NWeek) :: Rain     ! Meteo Input Data for Baseline
! SPI
Real, Parameter :: amssng = -99.0
!Real :: beta_spi(NMonth), gamm(NMonth), pzero(NMonth), SPI(Nst_Max,NYear*NWeek,NLen)
!=====================
Real :: beta_spi(NWeek), gamm(NWeek), pzero(NWeek), SPI(Nst_Max,NYear*NWeek,NLen)
!=====================
!Thong tin tram
Integer :: NSt, idx, Nlon, Nlat, Ntime, Nt_slice, idecade, Nmon_week, Nst_Fcs
Real, Dimension (Nst_Max) :: H, Lat, Lon !, AWC, SSMax, SUMax
Real :: DatBL(NYear*NWeek, Nst_Max) ! ivar=1: Tx, ivar=2: Tm, ivar=3: Rain
Integer :: ivar, ist, it, N, iy, j, it1, it2, iyr
Real :: pl, pmu, ts1, ms1, ts2, ms2, ps, zedry1, zewet1
!------------------------------------------------------------------------------
! Cho hang nam va ca chuoi thoi gian
Real, Dimension(Nst_Max,Nyear,NLen) :: NDYrs ! so thang han cua tung nam ung voi cac qui mo thoi gian
Real, Dimension(Nst_Max,Nyear,NLen) :: SumSPIYrs ! tong SPI<Nguong cua tung nam ung voi cac qui mo thoi gian
Real, Dimension(Nst_Max,Nyear,NLen) :: SPI_3m ! tong SPI cua 3 thang nao do trong nam (vi du 2-3-4)

Real, Dimension(Nyear,NLen) :: Area_Ext1_2,Area_Ext1_3,Area_Ext1_4,Area_Ext3_234  ! Tong so o luoi/tram co Han
!                              Thang 2, 3, 4 va Ba thang 2-3-4
! Cho ca giai doan:
Real :: DrEvent(Nst_Max,NYear*NWeek,NLen) ! chua gia tri 0/1 cua SPI<Nguong ung voi cac qui mo thoi gian
Real :: Spell_Dura(Nst_Max,0:NYear*NWeek,NLen) ! Chua so dot va do dai cac dot han
                                    ! Spell_Dura(0,NLen) la So dot han cua tung qui mo thoi gian
Real :: ND_All_Yrs(Nst_Max,NLen), SumSPI_All_Yrs(Nst_Max,NLen) ! So thang han vaf Tong SPI ca giai doan
! Sen's Slopes
Real, Dimension(Nst_Max,NLen) :: SenSl_NDYrs,Prob_NDYrs,SenSl_SumSPIYrs,Prob_SumSPIYrs
Real, Dimension(Nst_Max,NLen) :: SenSl_SPI_3m,Prob_SPI_3m
Integer :: NumYear ! Do dai chuoi thoi gian chap nhan duoc
! For Grads:
Integer :: Nx, Ny
Character (Len=50) :: Grads_File, Field, DrScale, GrFile, FieldName
Real :: Lon0, Lat0, DelX
!------------------------------------------------------------------------------
! Cac tham so cho chuong trinh du bao
Integer :: YrFcs1, YrFcs2, ilen
!Real :: ZpalmerFcs(Nst_Max,NLengFcs),VplmFcs(Nst_Max,NLengFcs)
Character (Len=255) :: Fname, BL_Dir, Res_Dir, FScale
!#####################
!! Lay thong tin vietnam
Integer :: nxlon, nxlat, ilon, jlat
Real, allocatable, dimension (:)   :: loninfo, latinfo, Xlon, Xlat
Real :: VN_Land(Nst_Max)   ! Cac o luoi tren dat lien VN va thong tin cac vung khi hau: B1=2, B2=3,..., N3=8
Real, allocatable, dimension (:,:) :: vietnam ! thong tin cac vung khi hau: B1=2, B2=3,..., N3=8
Logical :: status, Append
Real :: Wk(Nst_Max)
!#####################
!------------------------------------------------------------------------------
DatBL(:,:) = -99. !NYear*NWeek, NVar, Nst_Max)
SPI(:,:,:) = -99.; 
NDYrs(:,:,:) = -99.; SumSPIYrs(:,:,:) = -99.; DrEvent(:,:,:) = -99; Spell_Dura(:,:,:) = -99
!------------------------------------------------------------------------------

Decade = "Month"
!Decade  = "Week"
!BL_Dir  = "/work/projects/G09/OBS_VN/Monthly/TXT_R_T2m_4_Drought/"; 
BL_Dir  = "/home/tanpv/SEACLID/Data/"
Res_Dir = "../RES/SPI_Res/"

if (trim(Decade) == "Month") then
   Fname = "R_APH_Mon_1961-2006.txt"
   FScale = "_Mon_"
   Nmon_week = NMonth; Ntime = NYear*NMonth
else
   Fname = "R_Week_from_Day_4_Drought_1981-2014.txt"   ! Tuan 7 ngay
   FScale = "_Week_"
   Nmon_week = NWeek; Ntime = NYear*NWeek
endif

!Print*,"Reading data ..."
!Call StationInfo (Station, Lon, Lat, H, Nst)
!print*,"Nst=",Nst
!print*, (ist," ",station(ist),ist=1,Nst)
!pause

! 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")
j=0
Do jlat=1,Nxlat
   Do ilon=1,Nxlon
      j=j+1; VN_Land(j) = vietnam(ilon,jlat)
   Enddo
Enddo
!print*,"nxlon,nxlat=",nxlon,nxlat
!print*,"Lon=",loninfo
!print*,"Lat=",latinfo
Do jlat=nxlat, 1, -1
!   print "(40F4.1)",(VN_Land(ilon,jlat),ilon=1,nxlon)
Enddo
!pause

!Call Read_data_Baseline
Call Read_data_netcdf
print*,"Outputs from Read_data_netcdf Nst=",Nst
print "(A,70F6.1)","Lon=", Lon(1:Nst)
print "(A,70F6.1)","Lat=", Lat(1:Nst)

!print*, (ist," ",station(ist),ist=1,Nst)

pause

! SPI
Print*,"Calculating, Writing and Ploting SPI ..."
Call Cal_SPI
Call Cal_Drou_Characs
! Obtain here: spi(Nst,Ntime,NLen)  SPI(Nst_Max,NYear*NWeek,NLen)
Do ilen=1,NLen ! ist=1,Nst 
!   Do ist=1,Nst ! j=1,NLen
!      if (.not.Empty_Ar1 (spi(ist,1:Ntime,j), Ntime)) then
!          print "(12F6.1)", (DrEvent(ist,it,j), it=1,Ntime)
!          print*,Station(ist), lon(ist), lat(ist)
!          pause
!      endif
      print "(40F5.1)", (ND_All_Yrs(ist,ilen), ist=Nst,1,-1)
!   Enddo
!    pause
Enddo

!Call Cal_Drou_Characs
! Obtain here: ND_All_Yrs(ist,ilen), SumSPI_All_Yrs(ist,ilen)
!              Spell_Dura(ist,1:Ntime,ilen), Spell_Dura(ist,0,ilen)
!              NDYrs(ist,1:NYear,ilen), SumSPIYrs(ist,1:NYear,ilen), DrEvent(ist,1:Ntime,ilen)
! DrEvent(ist,1:Ntime,ilen)

Call Cal_Trend    ! Trend

print*,"Checking... Nlon, Nlat..", Nlon, Nlat
! Subroutine Write_Vector_Grid_Data_to_Grads (X,Nx,Ny,Grads_File, Field,DrScale, Lon0, Lat0,DelX)
Open (3, file="Drought_Time_series.txt",status="Unknown")
Call Write_Time_Series_to_TEXT (Area_Ext1_2, 3, "Dr_Ext_Mon_2")
Call Write_Time_Series_to_TEXT (Area_Ext1_3, 3, "Dr_Ext_Mon_3")
Call Write_Time_Series_to_TEXT (Area_Ext1_4, 3, "Dr_Ext_Mon_4")
Call Write_Time_Series_to_TEXT (Area_Ext3_234, 3, "Dr_Ext_Mon_234")

Grads_File = trim(Res_Dir)//"Drou_Charac.nc"; Append = .False. 
Do ilen=1, NLen 
   ! Number of Drought Months
   Write(Field,"(A,I2.2)") "N_mon_Dr_", Len_SPI(ilen);   Wk(1:Nst) = ND_All_Yrs(1:Nst,ilen)
   Call Write_2D_to_NC (Wk, trim(Grads_File), trim(Field), Append); Append = .True.

   ! Sum of SPI < Nguong
   Write(Field,"(A,I2.2)") "Sum_SPI_", Len_SPI(ilen);   Wk(1:Nst) = SumSPI_All_Yrs(1:Nst,ilen)
   Call Write_2D_to_NC (Wk, trim(Grads_File), trim(Field), Append); Append = .True.

   ! Spell_Dura(ist,0,ilen)
   Write(Field,"(A,I2.2)") "N_Spell_", Len_SPI(ilen);   Wk(1:Nst) = Spell_Dura(1:Nst,0,ilen)
   Call Write_2D_to_NC (Wk, trim(Grads_File), trim(Field), Append); Append = .True.

   ! Sen's Slope of Number of Drought months 
   Write(Field,"(A,I2.2)") "Sen_NDr_", Len_SPI(ilen);   Wk(1:Nst) = SenSl_NDYrs(1:Nst,ilen)
   Call Write_2D_to_NC (Wk, trim(Grads_File), trim(Field), Append); Append = .True.
   
   ! Sen's Slope of Sum SPI<Nguong in year
   Write(Field,"(A,I2.2)") "Sen_SPI_Yr", Len_SPI(ilen);   Wk(1:Nst) = SenSl_SumSPIYrs(1:Nst,ilen)
   Call Write_2D_to_NC (Wk, trim(Grads_File), trim(Field), Append); Append = .True.

   ! Sen's Slope of Sum SPI<Nguong in 3 months
   Write(Field,"(A,I2.2)") "Sen_SPI_3m", Len_SPI(ilen);   Wk(1:Nst) = SenSl_SPI_3m(1:Nst,ilen)
   Call Write_2D_to_NC (Wk, trim(Grads_File), trim(Field), Append); Append = .True.

   ! Geographical Extent:  Area_Ext(iyr,ilen)
!   Write(Field,"(A,I2.2)") "Sum_SPI_", Len_SPI(ilen);   Wk(1:Nst) = SumSPI_All_Yrs(1:Nst,ilen)
!   Call Write_2D_to_NC (Wk, trim(Grads_File), trim(Field), Append); Append = .True.

   ! Trend of No Dought months
! Real :: SenSl_NDYrs(Nst_Max,NLen),Prob_NDYrs(Nst_Max,NLen),SenSl_SumSPIYrs(Nst_Max,NLen),Prob_SumSPIYrs(Nst_Max,NLen)
!   Grads_File = "SenSl_NDYrs_"; Field = "SenNDYrs_"; write(DrScale,"(I2.2,A)") Len_SPI(ilen),"_Mon"
!   GrFile = trim(Grads_File)//trim(DrScale)//".nc"; FieldName = trim(Field)//trim(DrScale)
!   Call Write_Vector_Grid_Data_to_Grads &
!        (SenSl_NDYrs(1:Nst,ilen), Nx, Ny, trim(GrFile), trim(FieldName), trim(DrScale), Lon0, Lat0, DelX)
   ! Trend of Sum of Dought IDX
!   Grads_File = "SenSl_SumSPIYrs_"; Field = "SenSumSPIYrs_"; write(DrScale,"(I2.2,A)") Len_SPI(ilen),"_Mon"
!   GrFile = trim(Grads_File)//trim(DrScale)//".nc"; FieldName = trim(Field)//trim(DrScale)
!   Call Write_Vector_Grid_Data_to_Grads &
!        (SenSl_SumSPIYrs(1:Nst,ilen), Nx, Ny, trim(GrFile), trim(FieldName), trim(DrScale), Lon0, Lat0, DelX)
Enddo
Close (3)
!-----------------------------------------------------------------------
!! Call Plot_Hovm_SPI
!! Call Write_SPI_To_TEXT (trim(Res_Dir)//"SPI"//trim(FScale)//".txt")
! Calculate Drought Frequency
!! Print*,"Calculating Drought Frequency..."
!! Call Cal_Drought_Freq
!! Print*,"Ploting Frequency..."
!! Call Plot_Hovm_Freq (FScale)
!! Print*,"Writing Frequency..."
!! Call Write_Frequency_TEXT (Res_Dir, FScale)
Print*,"Done!"
Deallocate (loninfo,latinfo,vietnam)
!==============================================================================
CONTAINS

Subroutine Read_data_netcdf
Character (Len=255) :: FileName, VarName
!Real, Allocatable :: Xlon(:), Xlat(:), XTime(:), X(:,:,:)
Real, Allocatable :: XTime(:), X(:,:,:)

Integer :: Nlon1, Nlat1, Nt, i, j, ist, it
  FileName="/home/tanpv/SEACLID/Data/APH_mm_hist_allyear_50y.nc"
  VarName ="pr"
  Call rd_ncdf_3d (trim(FileName),Xlon,Xlat,XTime, Nlon1, Nlat1, Nt, X, trim(VarName))

  print*,"In Read_data_netcdf: Nlon, Nlat, Nt=", Nlon1, Nlat1, Nt
  print*,Xlon, Xlat
!  pause
!  do it=1,Nt
!     do j=Nlat1,1, -1
!        print "(40F4.1)", (X(i,j,it), i=1,Nlon1)
!     enddo
!     pause
!  Enddo

  Where (X < 0) X = -99; 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)
        Do it=1, Nt
           DatBL(it,ist) = X(i,j,it)
     !      print*,DatBL(it,ist)
        Enddo
     Enddo
  Enddo
  Nst = ist; Ntime = Nt; Nlon = Nlon1; Nlat = Nlat1
  print*,"Nst,Ntime",Nst,Ntime
  pause
End Subroutine Read_data_netcdf
!------------------------------------------------------------------------------
Subroutine Print_Test 
Integer :: it, ist
Do it=1, Ntime
   write (*,"(20F7.1)") (DatBL(it, ist),ist=1, Nst)
Enddo
End Subroutine Print_Test
!------------------------------------------------------------------------------
Subroutine Read_data_test ! Data from Tan
Integer :: iyt, ist, iyr, year, imw, mon, N
Character (Len=21) :: Stname
Real :: w(Nmon_week), xlon, xlat, HH
   print*, Nst, Ntime
   Open (1,file="/home/tanpv/PDSI/SPI/src_new/stdin")
   Do ist=1, Nst
      read (1, "(A)") Stname
      read (1,*) xlon, xlat, HH
      
      Do N=1, Ntime
         read (1, *) year, Mon, DatBL(N,ist)
         print*, year, Mon, DatBL(N,ist)
         DatBL(N,ist) = DatBL(N,ist)/100.
      Enddo
   Enddo
End Subroutine Read_data_test
!------------------------------------------------------------------------------

Subroutine Read_data_Baseline ! Data from Tan
Integer :: iyt, ist, iyr, year, imw, N, IO
Character (Len=21) :: Stname
Real :: w(Nmon_week)

   Open (1,file=trim(BL_Dir)//trim(Fname))
   ist = 1; IO = 0
   Do ! ist=1, Nst
      read (1, '(A)',IOStat=IO) Stname
      if (IO /= 0) exit
      Station(ist) = trim(Stname)
      read (1,"(2F8.3,F8.1)",IOStat=IO) lon(ist), lat(ist) !, H(ist)
      if (IO /= 0) exit
      N=0
      Do iyr=YBas1, YBas2
         read (1,*) year, (w(imw), imw=1,Nmon_week)
         Do imw=1, Nmon_week
            N=N+1; DatBL(N,ist) = w(imw)/100.
         Enddo
      Enddo
      ist = ist + 1
   Enddo
   Nst = ist - 1
End Subroutine Read_data_Baseline
!------------------------------------------------------------------------------
Real Function HSTQ (X, Y, N)
Integer :: N, i, j
Real :: X(N), Y(N), tbx, tby, sx, sy, rxy
tbx=0.; tby=0.; sx=0.; sy=0.;  rxy=0.
Do i=1, N
   tbx=tbx+x(i); tby=tby+y(i); sx=sx+x(i)*x(i); sy=sy+y(i)*y(i); rxy=rxy+x(i)*y(i)
Enddo
tbx=tbx/real(n); tby=tby/real(n); sx=sx/real(n); sy=sy/real(n); rxy=rxy/real(n)
sx=sqrt(sx-tbx*tbx); sy=sqrt(sy-tby*tby); rxy=rxy-tbx*tby 
HSTQ = rxy/(sx*sy)
End Function HSTQ
!------------------------------------------------------------------------------
Subroutine StationInfo (StaName,Lon, Lat, H, Nst) ! Reg, Nst)
Implicit None
Integer, Parameter :: NSta = 1000
Real :: Lon(NSta), Lat(NSta), H(NSta)
Integer :: i, i1, i2, i3, i4, Reg(NSta), IO, ist, Nst, L
Character (Len=21) StaName(NSta), st1, st2, st3, st4,st5
Character :: Str*255, Data_Dir*255

Data_Dir="./"
Open (1,file=Trim(Data_Dir)//"Danh_Sach_Tram_KT_N2.txt") ! Ten file danh sach tram day du
read (1,'(A)') Str   ! Bo dong dau
IO = 0; ist = 0
Do
   read (1,'(A)',IOstat = IO) Str
   if (io /= 0.or.trim(Str)=="") exit
   L = Len(trim(Str))
   st1=""; st2=""; st3=""; st4=""; st5=""
   i1=0; i2=0; i3=0; i4=0; i=0
   do while (i < L)  ! Danh dau so thu tu cua dau TAB trong xau
      i=i+1
      if (i1==0.and.Str(i:i)==Char(9)) then
         i1=i
      else if (i2==0.and.Str(i:i)==Char(9)) then
         i2=i
      else if (i3==0.and.Str(i:i)==Char(9)) then
         i3=i
      else if (i4==0.and.Str(i:i)==Char(9)) then
         i4=i
      endif
   enddo
   st1=Str(1:i1-1); st2=Str(i1+1:i2-1); st3=Str(i2+1:i3-1); st4=Str(i3+1:i4-1); st5=Str(i4+1:L)
   ist = ist+1  ! So thu tu trong danh sach tram
   StaName(ist) = st1   ! Tem tram
   read (st2,*) Lon(ist) ! Kinh do
   read (st3,*) Lat(ist) ! Vi do
   read (st4,*) H(ist)   ! Do cao
   read (st5,*) Reg(ist) ! Ky hieu vung
!   print '(I3,A,3F8.3,I3)',ist,Station(ist),Lon(ist), Lat(ist), H(ist), Reg(ist)
Enddo
Nst = ist  ! Tong so tram trong danh sach tram
Close (1)
End Subroutine StationInfo
!------------------------------------------------------------------------------
Logical Function Thang_Han (Mon,Rain)
! Xac dinh thang han cho Trung Bo (N.T.Hieu)
Integer :: Mon
Real :: Rain
Logical :: Han
Han = .false.
Select Case (Mon)
   Case (11, 12, 1, 2)
      if (Rain < 20) Han = .true.
   Case (9, 10, 3, 4)
      if (Rain < 30) Han = .true.
   Case (5, 6, 7, 8)
      if (Rain < 80) Han = .true.
End Select
Thang_Han = Han
End Function Thang_Han
!------------------------------------------------------------------------------
Subroutine Contingency_table (Rain, Dr_idx, N, Thres, FBI, PC, POD, HSS)
Integer :: N, i, Mon, Obs, Fcs
Real :: Rain(N), Dr_idx(N), Thres, a, b, c, d, FBI, PC, POD, HSS
a = 0; b = 0; c = 0; d = 0
Do i=1, N
   Mon = mod(i, 12); if (Mon == 0) Mon = 12
   if (Thang_Han(Mon,Rain(i))) then; Obs = 1; else; Obs = 0; endif
   if (Dr_idx(i) <= Thres) then; Fcs = 1; else; Fcs = 0; endif
   if (Obs==1.AND.Fcs==1) a = a + 1
   if (Obs==0.AND.Fcs==1) b = b + 1
   if (Obs==1.AND.Fcs==0) c = c + 1
   if (Obs==0.AND.Fcs==0) d = d + 1
Enddo
FBI=(a+b)/(a+c); PC=(a+d)/N; POD=a/(a+c); HSS=2*(a*d-b*c)/((a+c)*(c+d)+(a+d)*(b+d))
End Subroutine Contingency_table
!------------------------------------------------------------------------------
Include "Sub_SPI.f90" ! Chuong trinh tinh chi so han SPI
Include "Sub_Out_SPI.f90"
Include "Sub_Dr_Charac.f90"
Include "../Sub_Trend_MK_Linear.f90"
Include "sub_wr_ncdf.F90"
Include "sub_rd_ncdf.F90"
!==============================================================================
END   ! Program