Program PDSI_SEACLID 
use netcdf
Real, parameter :: amssng = -99.
Integer, Parameter :: NVar = 2, Nst_Max = 5000   ! So Grid Cells hoac so tram
Integer, Parameter :: NMonth = 12
Integer, Parameter :: YBas1 = 1966, YBas2 = 2005, NYear = YBas2-YBas1+1, maxyrs = NYear
Character (Len=4) :: YT(NVar) = (/"T2m","R"/)
Integer, Parameter :: NPDI = 3 ! 1=PDSI, 2=

Character (Len=21) :: Station(Nst_Max)
Real, Dimension (NVar, Nst_Max, YBas1:YBas2, NMonth) :: Dat_Mon ! ivar=1: T2m ivar=2: Rain
Real, Dimension (NVar, Nst_Max, NMonth) :: DatNor               ! Normal Climate Condition
Real :: AWC(Nst_Max), xlon(Nst_Max), xlat(Nst_Max)
! Outputs: 
Real, Dimension (Nst_Max,NYear*NMonth) :: mon_ori_pdsi, mon_ori_wplm, mon_ori_zind
Real, Dimension (Nst_Max,NYear*NMonth) :: mon_sel_pdsi, mon_sel_wplm, mon_sel_zind
Real, Dimension (Nst_Max,NYear*NMonth) :: Rain_Mon
!
Integer, Parameter :: Nidx=7, Ntq = 7   ! Kich thuoc Ma tran Tuong quan (1..3: Mon_Ori; 4..6: Mon_Sel; 7: Rain
Real, Dimension (Nst_Max,NYear*NMonth,Ntq) :: PMon  ! Luu tat ca cac chuoi PDSI de tinh HSTQ
Real, Dimension (Nst_Max,Ntq,Ntq) :: TQ_M                ! Ma tran HSTQ Mon & Week
!------------------------------------------------------------------------------
Character (Len=255) :: Fname, BL_Dir, Dat_Dir, Res_Dir
Character (Len=30) :: FMT
Logical :: status, Empty

Integer :: Nst, Nlon, Nlat, IO, Ntime, it
Integer :: ist, iyr, imw, ivar, isys, system
!! Lay thong tin vietnam
Integer :: nxlon, nxlat, ilon, jlat
Real, allocatable, dimension (:)   :: loninfo, latinfo, lon, lat
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, Xlonlat
!------------------------------------------------------------------------------
Dat_Mon(:,:,:,:) = -99.; DatNor(:,:,:) = 0.; Rain_Mon(:,:)=-99
!All_IDX(:,:,:) = -99; Leng_Spel(:,:,:,:) = -99; Drought(:,:,:,:) = -99
BL_Dir  = "/home/tanpv/SEACLID/Data/"
Res_Dir = "/home/tanpv/SEACLID/PDSI/my_result/"

Print*, "Reading Meteo data and station info ...."
FName="/home/tanpv/SEACLID/Data/APH_hist_allyear_40y.nc"  ! APH data 1966-2005
Call Read_data_netcdf (trim(FName))

Do ist=1, Nst
   it=0
   Do iyr=YBas1, YBas2
      Do imw=1, NMonth
         it=it+1
         Rain_Mon(ist,it) = Dat_Mon(2,ist,iyr,imw)
      Enddo
   Enddo
Enddo
Call putAWC(AWC, xlon, xlat, Nst); Where (AWC<0) AWC=6.7 ! Assign to gridboxes located along coast
!Call Print_test (AWC, .True.)  
!pause
!Call Print_test (vietnam, .False.)
!pause

isys = system ("rm -Rf /home/tanpv/SEACLID/PDSI/my_result")
mon_ori_pdsi(:,:)=-99; mon_ori_wplm(:,:)=-99; mon_ori_zind(:,:)=-99; 
mon_sel_pdsi(:,:)=-99; mon_sel_wplm(:,:)=-99; mon_sel_zind(:,:)=-99;

Do ist=1, Nst
   if (VN_Land(ist)>=1.AND.VN_Land(ist)<=7) then
  !    Empty = Empty_Stn (Dat_Mon, ist)
  !  if (Empty) then
  !      print*,"Khuyet SL: "//trim(station(ist))
  !  else
      Print*, "Calculation for station: "//trim(station(ist))
      Call Write_to_PDSI ("monthly", Dat_Mon,  NMonth, ist)

      isys = system ("./pdsi ")
      ! PostProcessing for station
      Call Read_Table (trim(Res_Dir)//"monthly/original/"//"PDSI.tbl", NMonth, mon_ori_pdsi(ist,:))
      Call Read_Table (trim(Res_Dir)//"monthly/original/"//"WPLM.tbl", NMonth, mon_ori_wplm(ist,:))
      Call Read_Table (trim(Res_Dir)//"monthly/original/"//"ZIND.tbl", NMonth, mon_ori_zind(ist,:))
      Call Read_Table (trim(Res_Dir)//"monthly/self_cal/"//"PDSI.tbl", NMonth, mon_sel_pdsi(ist,:))
      Call Read_Table (trim(Res_Dir)//"monthly/self_cal/"//"WPLM.tbl", NMonth, mon_sel_wplm(ist,:))
      Call Read_Table (trim(Res_Dir)//"monthly/self_cal/"//"ZIND.tbl", NMonth, mon_sel_zind(ist,:))

      PMon(ist,:,1)=mon_ori_pdsi(ist,:); PMon(ist,:,2)=mon_ori_wplm(ist,:); PMon(ist,:,3)=mon_ori_zind(ist,:)
      PMon(ist,:,4)=mon_sel_pdsi(ist,:); PMon(ist,:,5)=mon_sel_wplm(ist,:); PMon(ist,:,6)=mon_sel_zind(ist,:)
      PMon(ist,:,7)=Rain_Mon(ist,:)
      Call HSTQ (PMon, TQ_M, ist)

  !  endif
   endif
Enddo
print*,"Calculating PDSI: Done !"

Fname = "/home/tanpv/SEACLID/PDSI/Results/PDSI.nc"
Call Write_PDSI_to_NetCDF(Fname)

!Call Write_HSTQ ("HSTQ_Mon.txt", TQ_M)  ! Not necessary

CONTAINS
!------------------------------------------------------------------------------
Subroutine Write_HSTQ (Fname, X)
Real, Dimension (Nst_Max,Ntq,Ntq) :: X
Integer :: ist, j,k
Character (Len=*) :: Fname
Character (Len=2) :: Title(Ntq) = (/"op","ow","oz","sp","sw","sz","Pr"/)
Open (3, file=trim(Fname), Status="Unknown")
write (3, "(A)") "              "//Fname
Do ist=1, Nst
   write (3,"(A40)") trim(Station(ist))
   write (3,"(8A8)") " ",(Title(j), j=1,Ntq)
   Do j=1,Ntq
      write (3,"(A4,4x,7F8.4)") Title(j),(X(ist,j,k),k=1,Ntq)
   Enddo
Enddo
Close (3)
End Subroutine Write_HSTQ
!------------------------------------------------------------------------------
Subroutine HSTQ (X, Y, ist)
Real, Dimension (Nst_Max,NYear*NMonth,Ntq) :: X  ! Luu tat ca cac chuoi PDSI de tinh HSTQ
Real, Dimension (Nst_Max,Ntq,Ntq) :: Y
Integer :: ist, N, it, ic, j, k
Real :: TBj, TBk, Sj, Sk, Rjk, tmp
N = NYear*NMonth
Do j=1, Ntq
   Do k=1, Ntq
      TBj=0; TBk=0; Sj=0; Sk=0; Rjk=0; ic=0
      Do it=1, N
         if (X(ist,it,j) /= -99. .AND. X(ist,it,k) /= -99.) then
            ic=ic+1
            TBj=TBj+X(ist,it,j); TBk=TBk+X(ist,it,k); Sj=Sj+X(ist,it,j)*X(ist,it,j)
            Sk=Sk+X(ist,it,k)*X(ist,it,k); Rjk=Rjk + X(ist,it,j)*X(ist,it,k)
         endif
      Enddo
      tmp = ic
      TBj=TBj/tmp; TBk=TBk/tmp; Sj=Sj/tmp-TBj*TBj; Sk=Sk/tmp-TBk*TBk; Rjk=Rjk/tmp
      Y(ist,j,k) = (Rjk-TBj*TBk)/SQRT(Sj*Sk)
   Enddo
Enddo
End Subroutine HSTQ
!------------------------------------------------------------------------------
Subroutine Read_Table (Fname, Nmon_week, Outp) 
Integer :: Nmon_week, it, iyr, imw, year, isys
Real :: wk(Nmon_week), Outp(Nmon_week*NYear)
Character (Len=*) :: Fname
Logical :: exists
Inquire (File=trim(Fname),exist=exists)
if (.Not.exists) return
Open (1, File=trim(Fname), status="Old")
it = 0
!print*,trim(Fname), "<<<<<<================"
Do iyr=YBas1, YBas2
   read (1, *) year, (wk(imw), imw=1, Nmon_week)
!   print*, year, (wk(imw), imw=1, Nmon_week)
   Do imw=1, Nmon_week
      it = it+1; Outp(it) = wk(imw)
   Enddo
Enddo
isys = system ("rm -f "//trim(Fname))
End Subroutine Read_Table
!------------------------------------------------------------------------------
Subroutine Write_to_PDSI (Dat_Dir, DatBL, Nmon_week, ist)
Integer :: ivar, ist, iyr, imw, isys, system
Real :: Tmp
Character (Len=*) :: Dat_Dir
Character (Len=255) :: Fname
Real, Dimension (NVar, Nst_Max, YBas1:YBas2, NMonth) :: DatBL
! Write data of one station to calculate PDSI
   Fname = "/home/tanpv/SEACLID/PDSI/my_dat"
   isys = system ('mkdir -p '//trim(Fname))
   Open (1, file=trim(Fname)//"/"//"parameter", status="Unknown")
   write (1,*) AWC(ist), xlat(ist)
   Do ivar=1, NVar
      if (ivar == 1) then
         Open (3, file=trim(Fname)//"/"//trim(Dat_Dir)//"_T", status="Unknown")
         Open (4, file=trim(Fname)//"/"//Dat_Dir(1:3)//"_T_normal",status="Unknown")
      else
         Open (3, file=trim(Fname)//"/"//trim(Dat_Dir)//"_P", status="Unknown")
         Open (4, file=trim(Fname)//"/"//Dat_Dir(1:3)//"_P_normal",status="Unknown")
      endif
      Do imw=1, Nmon_week   ! Tinh trung binh khi hau
         Tmp = 0
         Do iyr=YBas1, YBas2
            Tmp = Tmp + DatBL(ivar,ist,iyr,imw)
         Enddo
         DatNor(ivar,ist,imw) = Tmp/Real(NYear)
      Enddo
      write (4,"(52F8.3)") (DatNor(ivar,ist,imw), imw=1, Nmon_week)
      Do iyr=YBas1, YBas2
         write (3,"(I4,52F8.3)") iyr,(DatBL(ivar,ist,iyr,imw), imw=1, Nmon_week)
!         write (*,"(I4,52F8.3)") iyr,(DatBL(ivar,ist,iyr,imw), imw=1, Nmon_week)
      Enddo
      Close (3); Close (4)
   Enddo
   Close (1)
End Subroutine Write_to_PDSI
!------------------------------------------------------------------------------
Subroutine Read_data_netcdf (FileName)
Character (Len=*) :: FileName
Real, Allocatable :: XTime(:), T2m(:,:,:), Pr(:,:,:) !, VN(:,:)
Real, allocatable, dimension (:) :: loninfo, latinfo
Integer :: nxlon, nxlat, ilon, jlat
Integer :: Nlon1, Nlat1, Nt, i, j, ist, it, iyr,imw

!  FileName="/home/tanpv/SEACLID/Data/APH_mm_hist_allyear.nc"  ! APH data (1986-2005)
  Call rd_ncdf_3d (trim(FileName), lon, lat, XTime, Nlon1, Nlat1, Nt, Pr, "pr")  ! Rainfall
  Deallocate (lon, lat, XTime) 
  Call rd_ncdf_3d (trim(FileName), lon, lat, XTime, Nlon1, Nlat1, Nt, T2m, "ts") ! Temperature
  
  FileName="/home/tanpv/SEACLID/MASK/VN.grid_0.25.nc"
  Call rd_ncdf_2d (trim(FileName),loninfo,latinfo,nxlon,nxlat, vietnam, "all")
!  Allocate (vietnam(nxlon,nxlat)); vietnam(:,:) = VN(:,:)
  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
!  Call Print_test (VN, .False.); pause

  print*,"In Read_data_netcdf: Nlon, Nlat, Nt=", Nlon1, Nlat1, Nt
  print*,lon, lat

  Where (Pr < 0) Pr = -99; Where (T2m < 0) T2m = -99; ist = 0
  Do j=1, Nlat1
     Do i=1, Nlon1
        ist=ist+1; write(Station(ist),"(A,I4.4)") "Point-",ist
        XLon(ist)=lon(i); XLat(ist)=lat(j);
        VN_Land(ist)=vietnam(i,j)
        it=0
        Do iyr=YBas1, YBas2  ! 1, NYear
           Do imw=1, Nmonth
              it=it+1
              if (T2m(i,j,it)/=-99) Dat_Mon(1,ist,iyr,imw) = T2m(i,j,it)*1.8 + 32.    ! Doi sang do F
              if (Pr (i,j,it)/=-99) Dat_Mon(2,ist,iyr,imw) = Pr (i,j,it)/25.4         ! Doi sang Inch
           Enddo
        Enddo
     Enddo
  Enddo
  Deallocate ( XTime, loninfo,latinfo) 
  Nst = ist; Ntime = Nt; Nlon = Nlon1; Nlat = Nlat1
  print*,"Nst,Ntime",Nst,Ntime
End Subroutine Read_data_netcdf
!------------------------------------------------------------------------------
Subroutine Write_PDSI_to_NetCDF(Fname)
Integer :: ist, jlat, ilon, it, ilen, ix
Real, Dimension (Nlon, Nlat, Ntime) :: X1, X2, X3, X4, X5, X6
Real, Dimension (Ntime) :: t
Character (Len=30) :: t_unit
Character (Len=*) :: Fname

Do it=1, Ntime
   t(it) = it-1
Enddo

ist = 0
Do jlat=1, Nlat
   Do ilon=1, Nlon
      ist = ist +1
      Do it = 1, Ntime
         X1(ilon,jlat,it) = mon_ori_pdsi(ist,it)
         X2(ilon,jlat,it) = mon_ori_wplm(ist,it)
         X3(ilon,jlat,it) = mon_ori_zind(ist,it)
         X4(ilon,jlat,it) = mon_sel_pdsi(ist,it)
         X5(ilon,jlat,it) = mon_sel_wplm(ist,it)
         X6(ilon,jlat,it) = mon_sel_zind(ist,it)
      Enddo
   Enddo
Enddo

t_unit = "months since 1986-01-15"
Call wr_ncdf_3d (trim(Fname), Lon, Lat, t, Nlon, Nlat, Ntime, X1, "ori_pdsi", trim(t_unit), amssng)
Call wr_ncdf_3d_apd (trim(Fname), X2, "ori_wplm", Nlon, Nlat, Ntime, amssng)
Call wr_ncdf_3d_apd (trim(Fname), X3, "ori_zind", Nlon, Nlat, Ntime, amssng)
Call wr_ncdf_3d_apd (trim(Fname), X4, "sel_pdsi", Nlon, Nlat, Ntime, amssng)
Call wr_ncdf_3d_apd (trim(Fname), X5, "sel_wplm", Nlon, Nlat, Ntime, amssng)
Call wr_ncdf_3d_apd (trim(Fname), X6, "sel_zind", Nlon, Nlat, Ntime, amssng)
Call wr_ncdf_2d_apd (trim(Fname), vietnam, "Subreg", Nlon, Nlat, amssng)

End Subroutine Write_PDSI_to_NetCDF
!------------------------------------------------------------------------------
SUBROUTINE putAWC(dat,lon,lat,nst)
   integer :: nst, ist
   real :: dat(nst), lon(nst), lat(nst)
   integer, parameter :: ncols = 1099, nrows = 549
   integer :: icol, irow, itmp
   character (len = *), parameter :: fileAWC = "/work/users/thanhnx/Func_Sub/AWC_Global.asc"
   real, parameter :: xcor = -180.0, ycor = -89.9056530001
   real, parameter :: cellsiz = 0.327697, iflag = -9999, oflag = -99
   real :: indat(ncols,nrows), x(ncols), y(nrows)
! <-~~~=|==-----------------------------------------------------------<o
! <-~~~=|== Read asc AWC file
   open (1,file = fileAWC,status = "old")
      ! <-~~~=| Skip reading 7 lines
      read (1,*) ; read (1,*) ; read (1,*) ; read (1,*)
      read (1,*) ; read (1,*) ; read (1,*)
      ! <-~~~=|
      do irow = 1, nrows
         itmp = nrows - irow + 1
         read (1,*) (indat(icol,itmp),icol = 1, ncols)
      enddo
   close (1)
   do icol = 1, ncols
      x(icol) = xcor + (icol - 1) * cellsiz
   enddo
   do irow = 1, nrows
      y(irow) = ycor + (irow - 1) * cellsiz
   enddo
   ! <-~~~=| Compute for new grid
   do ist = 1, nst
      icol = nrtpos(x,lon(ist),ncols,oflag)
      irow = nrtpos(y,lat(ist),nrows,oflag)
      dat(ist) = indat(icol,irow)*0.0393701
      if (dat(ist) == iflag) then
         dat(ist) = oflag   ! Doi -99.9 --> 99.0 
      endif
   enddo
return
END SUBROUTINE putAWC
!--------------------------------------------------------------------
integer FUNCTION nrtpos(x,y,n,flag)
   integer :: n, i, itmp
   real :: flag, x(n), y, minx

   minx = 999999
   itmp = 0
   do i = 1, n
      if (x(i) /= flag .and. y /= flag) then
         if (abs(x(i) - y) < minx) then
            minx = abs(x(i) - y)
            itmp = i
         endif
      endif
   enddo

   nrtpos = itmp
   return
END FUNCTION nrtpos
!--------------------------------------------------------------------
Subroutine Print_test (X, One_D)  ! Print test of 1-D array X(Nst)
Integer :: ist, jlat, ilon
Real :: X(Nst_Max)
Logical :: One_D
!if (One_D) then
!   ist=0; if (.not.Allocated (Xlonlat)) Allocate(Xlonlat(Nlon,Nlat))
!   Do jlat=1,Nlat
!      Do ilon=1,Nlon
!         ist=ist+1; Xlonlat(ilon,jlat)=X(ist) ! Rain_Mon(ist,1)
!      enddo
!   Enddo
!endif
Do jlat=Nlat,1,-1
   print "(40F6.1)", (Xlonlat(ilon,jlat),ilon=1,Nlon)
Enddo
!pause ; 
!if (One_D) Deallocate(Xlonlat)
End Subroutine Print_test
!------------------------------------------------------------------------------
Logical Function Empty_Stn (DatBL, ist)
! Tra ve .True. neu tram ist khuyet SL
Real, Dimension (NVar, Nst_Max, YBas1:YBas2, NMonth) :: DatBL
Integer :: N, ist,iyr,imon, k1, k2
Logical :: L1, L2

L1=.False.; N=0
Do iyr=YBas1, YBas2
   Do imon=1,NMonth
      if (DatBL(1,ist,iyr,imon)<0.OR.DatBL(2,ist,iyr,imon)<0) then
         N = N+1
         if (N>10) then ; L1=.True.; Exit; endif 
      endif
   Enddo
Enddo
Empty_Stn = L1 !.OR.L2    ! T hoac P khuyet SL
End Function Empty_Stn
!--------------------------------------------------------------
Include "sub_rd_ncdf.F90"
Include "sub_wr_ncdf.F90"
!Include "Sub_Plot.f90"
!Include "Sen_Slope.f90"
!--------------------------------------------------------------------
END ! Program