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

Open (1, file = "File_List.txt"); Read (1, Filenames)

Do ifile=1, Nfile
   Dat_Mon(:,:,:,:) = -99.; DatNor(:,:,:) = 0.; Rain_Mon(:,:)=-99
   Res_Dir = "/home/tanpv/SEACLID_New/PDSI/my_result/"

   Print*, "Reading Meteo data and station info ...."
   FName=Fname_inp(ifile)    ! Data file to be read
   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_New/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 !"
   print*,"Writing Results to file:", Fname_Out(ifile)
   Fname = Fname_Out(ifile) !"/home/tanpv/SEACLID_New/PDSI/Results/PDSI_APH_BL.nc"
   Call Write_PDSI_to_NetCDF(Fname)

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

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_New/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

  if (Allocated(lon)) Deallocate(lon); if (Allocated(lat)) Deallocate(lat)
  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, "tas") ! Temperature
  
  FileName="/home/tanpv/SEACLID_New/MASK/VN.grid_0.25.nc"
  if (Allocated(vietnam)) Deallocate(vietnam)
  Call rd_ncdf_2d (trim(FileName),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
!  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)*DayOfMon(imw,iyr)/25.4         ! Doi sang Monthly Total in Inch
!print*, iyr,imw,it
           Enddo
        Enddo
     Enddo
  Enddo
  Deallocate ( XTime, loninfo,latinfo) 
  Nst = ist; Ntime = Nt !-240; 
  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
!------------------------------------------------------------------------------
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_rd_ncdf.F90"
Include "sub_wr_ncdf.F90"
!Include "Sub_Plot.f90"
!Include "Sen_Slope.f90"
!--------------------------------------------------------------------
END ! Program