Program Dr_IDXs_SEACLID
  ! Calculating PDSI, SPI, SPEI 
  use netcdf

!------------------------------------------------------------------------------

  Real, parameter    :: amssng = -99., Miss=amssng
  Integer, Parameter :: NMonth = 12, Nmon=NMonth
  Integer, Parameter :: NVar = 2, Nreg = 7, Nmem = 7, Nlead = 6
  Integer, Parameter :: Ncase = 1       !  obs,mod,cli,ann
  Character (Len=3), Parameter :: Cases(Ncase)=(/"obs"/) !,"mod","cli","ann"/)

  Integer, Allocatable :: Len_SPI(:)
  !Real, Allocatable    :: DatBL(:,:), DatTMP(:,:), Rain(:) ! For calculating SPI, SPEI
  Real, Allocatable    :: Rain(:), T2m(:),xspi(:)

  Real, Allocatable    :: beta_spi(:), gamm(:), pzero(:), SPI(:,:,:), SPE(:,:,:)
  Real                 :: beta_spei(0:2), logLogisticParams(0:Nmon,0:2), DisPars(0:Nmon,0:2)

  Real, Allocatable :: Regs(:,:), Lon(:), Lat(:), xlon(:), xlat(:) 
  Real, Allocatable :: Pr_All(:,:,:,:,:,:), T2_All(:,:,:,:,:,:)   ! Nlon,Nlat,Nmem,Nlead,Ntime,Ncase
  Real, Allocatable :: xmem(:), xlead(:), xtime(:)
  Real, Allocatable :: SPEI_All(:,:,:,:)                  ! Nlon,Nlat,Ntime,Nlen
  Real, Allocatable :: SPI_All (:,:,:,:) 

  ! Notes:
  ! 1) Dat_Mon(ivar,ist,iyr,imw), where ivar=1:Temp (F), ivar=2:Rainfall (inch); ist=1..Nst; iyr=YBas1..YBas2; imw=1..Nmon (or Nweek) 
  !    Dat_Mon used to calculating PDSI
  ! 2) Temp_Inp (Temp in C deg), Rain_Inp (Rainfall in mm/month) (1:Nlon1,1:Nlat1,1:Ntime), Ntime=(YBas2-YBas1+1)*Nmon
  !    Temp_Inp and Rain_Inp to be used to write out as input data
  ! 3) Prec (mm/month), Temp (C) (1:Ntime,1:Nst)
  ! 4) DatBL (mm/month), DatTMP (C) (1:Ntime,1:Nst): Used to calculating SPI, SPEI

!------------------------------------------------------------------------------
  Character (Len=255) :: Fname, BL_Dir, R_Fname, T_Fname, Out_Dr_Fname
  Character (Len=255) :: Inp_Dir, Out_Dir, Mask_File
  Character (Len=3)   :: TVar, RVar
  Character (Len=10),Allocatable :: Model_nm(:)
  Logical :: status, Empty, use_baseline

  Integer :: Yr1, Yr2, Nmodel, NYear, maxyrs
  Integer :: NLen, Nst, Nlon, Nlat, IO, Ntime, Nt, it, Landpoint
  Integer :: imod, ist, iyr, imw, ivar, isys, system, ilon, ilat
  Integer :: Nt1
  Real    :: latx

!------------------------------------------------------------------------------
! For writing to NC file:
  Real, allocatable, Dimension (:,:,:,:) :: X4d
  Real, allocatable, Dimension (:) :: X1d,Xt
!------------------------------------------------------------------------------

  Namelist /File_Dir_Names/ Inp_Dir, Model_nm, Mask_File, Out_Dir, TVar, RVar
  Namelist /Data_Info/      Yr1, Yr2, Nmodel, NLen
  Namelist /Model_name/     Model_nm
  Namelist /Dr_Timescale/   Len_SPI
!------------------------------------------------------------------------------

  Open (1, file="Namelist_Dr_Reg.txt", Status="Old")     ! Tinh voi APH 1961-2015

  Read (1, File_Dir_Names)
  Read (1, Data_Info)
  Allocate (Len_SPI(NLen)); Allocate (Model_nm(Nmodel))
  Read (1, Model_name)
  Read (1, Dr_Timescale)
  Close (1)

  Yr1=1982; Yr2=2021     ! Tinh cho SL OBS
  Nyear = (Yr2 - Yr1 + 1); Ntime = Nyear*NMonth; maxyrs = NYear
  Nt = Ntime ! - Nmon   ! Neu co -Nmon thi tru SL nam cuoi cung (vi chi co cac thang du bao ung voi cac LT)
  Call Read_Regs 
  print*,"Parameters: Nlat,Nlon,Nst,Nlen,Ntime=",Nlat,Nlon,Nst,Nlen,Ntime
  imod=1
     Call Reinit
     Call Read_all_data (Model_nm(imod))
     print*,"Reading data done ..."
     !------------------------------------------------------------------------------
     print*,"Beginning: Ntime, Nst=",Ntime, Nst
     print*,"Parameters: Nlat,Nlon,Nst,Nlen,Ntime=",Nlat,Nlon,Nst,Nlen,Ntime

      Ntime1=Ntime           !!!!!!!!!!
      ilead=1; imem=1; icas=1       ! Do voi OBS tat ca cac model, cac mem va cac lead deu nhu nhau ==> Chi lay 1 tr/h
      use_baseline=.false.
      SPEI_All=Miss; SPI_All=Miss

      Do ilon=1,Nlon
         Do ilat=1,Nlat
            if (Regs(ilon,ilat)>0) then
               Rain(1:Ntime)=Pr_All(ilon,ilat,imem,ilead,1:Ntime,icas)
               T2m (1:Ntime)=T2_All(ilon,ilat,imem,ilead,1:Ntime,icas)
               latx=xlat(ilat)
               Nt = Min(Length_series (Rain,Ntime,Miss), Length_series (T2m,Ntime,Miss))
               Do ilen=1,Nlen
                  Call Calc_SPEI_One_Stn (T2m, Rain, latx, Nt, Nmon, Len_SPI(ilen), logLogisticParams, use_baseline, xspi, Miss)                  
                  SPEI_All(ilon,ilat,1:Nt,ilen)=xspi(1:Nt)
                  Call Cal_SPI_One_Stn (Rain,Nt,Len_SPI(ilen),alpha,beta_spi,gamm,pzero,xspi,use_baseline,Miss)
                  SPI_All (ilon,ilat,1:Nt,ilen)=xspi(1:Nt)
               Enddo
            endif
         Enddo
      Enddo 
      write (Fname,"(A,I4.4,A,I4.4,A)") trim(Out_Dir)//"SPEI_SPI_obs_",Yr1,"-",Yr2,".nc"
      Call Write_Dr_IDXs_to_file (Fname,ilead,Model_nm(imod)) 
   
   CONTAINS
   !------------------------------------------------------------------------------
   Subroutine Reinit
      if (Allocated(T2_All)) Deallocate (T2_All)
      if (Allocated(Pr_All)) Deallocate (Pr_All)

      if (Allocated (Rain))     Deallocate (Rain)
      if (Allocated (T2m))      Deallocate (T2m)
      if (Allocated (xspi))     Deallocate (xspi)
      !if (Allocated (DatBL))    Deallocate (DatBL)
      !if (Allocated (DatTMP))   Deallocate (DatTMP)
      if (Allocated (beta_spi)) Deallocate (beta_spi)
      if (Allocated (gamm))     Deallocate (gamm)
      if (Allocated (pzero))    Deallocate (pzero)
      if (Allocated (SPEI_All)) Deallocate (SPEI_All)
      if (Allocated (Lat))      Deallocate (Lat)
      if (Allocated (X1d))      Deallocate (X1d)

      !Allocate (Rain(Ntime),DatBL(Ntime,Nst),DatTMP(Ntime,Nst))
      Allocate (Rain(Ntime), T2m(Ntime), xspi(Ntime))
      Allocate (beta_spi(Nmonth), gamm(Nmonth), pzero(Nmonth))
      !Allocate (SPI(Nst,Ntime,NLen), SPE(Nst,Ntime,NLen))
      Allocate (SPEI_All(Nlon,Nlat,Ntime,Nlen))
      Allocate (SPI_All (Nlon,Nlat,Ntime,Nlen)) 
      Allocate (X1d(Ntime))
   End Subroutine Reinit
   !------------------------------------------------------------------------------
   Subroutine Write_Dr_IDXs_to_file (Fname,ilead,Model_nm)
     Character (Len=*) :: Fname,Model_nm
     Character (Len=250) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm
     Integer :: it, imem, ilead, ireg, icas, ilen, imod
     Real :: xlen(Nlen), tmp,tmp1, thres, xxx(Ntime)
     Real, Allocatable :: X41(:,:,:,:),X42(:,:,:,:),Y4d(:,:,:,:)

     t_unit="months since 1982-01-01 "; var_unit="-";
     lon_nm="lon"; lat_nm="lat"; z_nm="timescale"; t_nm="time"
     Do ilen=1,Nlen; xlen(ilen)=Len_SPI(ilen);  Enddo

     Allocate (X4d(Nlon,Nlat,Nlen,Ntime),Y4d(Nlon,Nlat,Nlen,Ntime))

     Do ilen=1,Nlen 
        Do ilon=1,Nlon
           Do ilat=1,Nlat
              Do it=1,Ntime
                 X4d(ilon,ilat,ilen,it)=SPEI_All(ilon,ilat,it,ilen)
                 Y4d(ilon,ilat,ilen,it)=SPI_All (ilon,ilat,it,ilen)
              Enddo
           Enddo
        Enddo
     Enddo  
     write(Var_ln,"(A,I2.2,A,I1)") "SPEI for observation, timescale = 1,3,6 months corresponding to z=1,2,3" 
     Call wr_ncdf_4d (Fname, xlon,xlat,xlen,xtime, Nlon,Nlat,Nlen,Ntime,X4d,"spei", &
                              t_unit,Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss)
     write(Var_ln,"(A,I2.2,A,I1)") "SPI for observation, timescale = 1,3,6 months corresponding to z=1,2,3"
     Call wr_ncdf_4d_apd (Fname,Y4d,"spi",Nlon,Nlat,Nlen,Ntime, &
                             Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss)
        !========================================
        ! Checking: No_missing_Ar1 (X,N,Miss)
        !Do ilen=1,Nlen
        !   ii=0; j=Len_SPI(ilen)-1
        !   Do ilon=1,Nlon
        !      Do ilat=1,Nlat
        !         if (Regs(ilon,ilat)>0) then
        !            xxx(1:Ntime-j)=X4d(ilon,ilat,ilen,j+1:Ntime)
        !            if (.not.No_missing_Ar1(xxx,Ntime-j,Miss)) then
        !               ii=ii+1; print*,ii
        !               print"(20F7.2)",xxx(1:Ntime-j)
        !            endif
        !         endif
        !      enddo
        !   Enddo
        !   pause
        !Enddo
        !========================================
        ! Xac dinh cac drought events:

        Allocate (X41(Nlon,Nlat,Nlen,Ntime),X42(Nlon,Nlat,Nlen,Ntime));

        thres=-0.5; X41=Miss
        Do ilon=1,Nlon
           Do ilat=1,Nlat
              Do ilen=1,Nlen
                 Do it=1,Ntime
                    tmp=X4d(ilon,ilat,ilen,it); tmp1=Y4d(ilon,ilat,ilen,it)
                    if (tmp/= Miss.AND.tmp <=thres) X41(ilon,ilat,ilen,it)=1
                    if (tmp/= Miss.AND.tmp >thres)  X41(ilon,ilat,ilen,it)=0
                    if (tmp1/=Miss.AND.tmp1<=thres) X42(ilon,ilat,ilen,it)=1
                    if (tmp1/=Miss.AND.tmp1>thres)  X42(ilon,ilat,ilen,it)=0
                 Enddo
              Enddo
           Enddo
        Enddo
        Var_ln="Observed mild drought events (dridx<=-0.5)"
        Call wr_ncdf_4d_apd (Fname,X41,"spei_mil",Nlon,Nlat,Nlen,Ntime, &
                             Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss)
        Call wr_ncdf_4d_apd (Fname,X42,"spi_mil",Nlon,Nlat,Nlen,Ntime, &
                             Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss)
        thres=-1.0; X41=Miss
        Do ilon=1,Nlon
           Do ilat=1,Nlat
              Do ilen=1,Nlen
                 Do it=1,Ntime
                    tmp=X4d(ilon,ilat,ilen,it); tmp1=Y4d(ilon,ilat,ilen,it)
                    if (tmp/=Miss.AND.tmp<=thres) X41(ilon,ilat,ilen,it)=1
                    if (tmp/=Miss.AND.tmp>thres)  X41(ilon,ilat,ilen,it)=0
                    if (tmp1/=Miss.AND.tmp1<=thres) X42(ilon,ilat,ilen,it)=1
                    if (tmp1/=Miss.AND.tmp1>thres)  X42(ilon,ilat,ilen,it)=0
                 Enddo
              Enddo
           Enddo
        Enddo
        Var_ln="Observed moderate drought events (dridx<=-1.0)"
        Call wr_ncdf_4d_apd (Fname,X41,"spei_mod",Nlon,Nlat,Nlen,Ntime, &
                             Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss)
        Call wr_ncdf_4d_apd (Fname,X42,"spi_mod",Nlon,Nlat,Nlen,Ntime, &
                             Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss)
        thres=-1.5; X41=Miss
        Do ilon=1,Nlon
           Do ilat=1,Nlat
              Do ilen=1,Nlen
                 Do it=1,Ntime
                    tmp=X4d(ilon,ilat,ilen,it); tmp1=Y4d(ilon,ilat,ilen,it)
                    if (tmp/=Miss.AND.tmp<=thres) X41(ilon,ilat,ilen,it)=1
                    if (tmp/=Miss.AND.tmp>thres)  X41(ilon,ilat,ilen,it)=0
                    if (tmp1/=Miss.AND.tmp1<=thres) X42(ilon,ilat,ilen,it)=1
                    if (tmp1/=Miss.AND.tmp1>thres)  X42(ilon,ilat,ilen,it)=0
                 Enddo
              Enddo
           Enddo
        Enddo
        Var_ln="Observed severe drought events (dridx<=-1.5)"
        Call wr_ncdf_4d_apd (Fname,X41,"spei_sev",Nlon,Nlat,Nlen,Ntime, &
                             Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss)
        Call wr_ncdf_4d_apd (Fname,X42,"spi_sev",Nlon,Nlat,Nlen,Ntime, &
                             Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss)
        thres=-2.0; X41=Miss
        Do ilon=1,Nlon
           Do ilat=1,Nlat
              Do ilen=1,Nlen
                 Do it=1,Ntime
                    tmp=X4d(ilon,ilat,ilen,it); tmp1=Y4d(ilon,ilat,ilen,it)
                    if (tmp/=Miss.AND.tmp<=thres) X41(ilon,ilat,ilen,it)=1
                    if (tmp/=Miss.AND.tmp>thres)  X41(ilon,ilat,ilen,it)=0
                    if (tmp1/=Miss.AND.tmp1<=thres) X42(ilon,ilat,ilen,it)=1
                    if (tmp1/=Miss.AND.tmp1>thres)  X42(ilon,ilat,ilen,it)=0
                 Enddo
              Enddo
           Enddo
        Enddo
        Var_ln="Observed extreme drought events (dridx<=-2.0)"
        Call wr_ncdf_4d_apd (Fname,X41,"spei_ext",Nlon,Nlat,Nlen,Ntime, &
                             Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss)
        Call wr_ncdf_4d_apd (Fname,X42,"spi_ext",Nlon,Nlat,Nlen,Ntime, &
                             Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss)
           !else
           !   Call wr_ncdf_4d_apd (Fname,X4d,Var_nm,Nlon,Nlat,Nmem,Ntime, &
           !                   Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss)
           !endif
     !   Enddo
     !Enddo
     Deallocate (X4d,X41,X42,Y4d)
   End Subroutine Write_Dr_IDXs_to_file
   !------------------------------------------------------------------------------
   Integer Function Length_series (X,Ninp,Miss)
     Integer :: Ninp,i,j
     Real :: X(Ninp), Miss

     j=0
     Do i=Ninp,1,-1
        if (X(i)==Miss) then
           j=j+1
        else
           exit
        endif
     Enddo
     !print*,Ninp,j
     Length_series=Ninp-j
   End Function Length_series
   !------------------------------------------------------------------------------
   Subroutine Read_Regs 
     Integer :: ilon,ilat, ireg, ic,j
     Real :: tmp

     Call rd_ncdf_2d (Mask_File,xlon,xlat,Nlon,Nlat,Regs,"reg")
     Nst=0;
     Do ilon=1,Nlon
        Do ilat=1,Nlat
           if (Int(Regs(ilon,ilat))>0) then
              Nst=Nst+1
           endif
        Enddo
     Enddo
     !print*,"Nst=",Nst; pause
   End Subroutine Read_Regs
   !------------------------------------------------------------------------------
   Subroutine Read_all_data (Model_nm)
      Character (Len=*) :: Model_nm
      Integer :: ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,icas, it1, it2
      Character (Len=255) :: Fname
      Real, Allocatable :: X4d(:,:,:,:), time(:), xlon(:),xlat(:),xmem(:) !,xtime(:)

      it1=Nmon+1      ! Tinh tu   1/1983
      it2=it1+Ntime-1 ! Tinh den 12/2021
      if (Allocated(X4d))  Deallocate (X4d);  if (Allocated(xmem))  Deallocate (xmem)
      if (Allocated(xlon)) Deallocate (xlon); if (Allocated(xlat))  Deallocate (xlat)
      if (Allocated(time)) Deallocate (time); if (Allocated(xtime)) Deallocate (xtime)
      Do ilead=1,Nlead
         write (Fname,"(A,I2.2,A)") trim(Inp_Dir)//trim(TVar)//"_"//trim(Model_nm)//"_lead_",ilead,".nc"
         print*,"Reading ",trim(Fname)
         Do icas=1,Ncase
            Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,trim(Cases(icas)))
            print*,"Nlo,Nla,Nmem,Ntime=",Nlo,Nla,Nme,Nt
            if (icas==1.AND.ilead==1) then 
               Allocate (T2_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase),time(Ntime)); 
               T2_All=Miss; time(1:Ntime)=xtime(1:Ntime)    ! Luu lai de write out
            endif
            Do imem=1,Nmem
               Do it=1,Ntime
                  T2_All(:,:,imem,ilead,it,icas)=X4d(:,:,imem,Nmon+it)  ! Luu vao T2_All voi time tu Nmon+1 = it1 (Tinh tu   1/1983)
               Enddo
            Enddo
            Deallocate (xlon,xlat,xmem,xtime,X4d)
         Enddo
         write (Fname,"(A,I2.2,A)") trim(Inp_Dir)//trim(RVar)//"_"//trim(Model_nm)//"_lead_",ilead,".nc"
         Do icas=1,Ncase
            Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,trim(Cases(icas)))
            print*,"Nlo,Nla,Nmem,Ntime=",Nlo,Nla,Nme,Nt
            if (icas==1.AND.ilead==1) then; Allocate (Pr_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase)); Pr_All=Miss; endif
            Do imem=1,Nmem
               Do it=1,Ntime
                  Pr_All(:,:,imem,ilead,it,icas)=X4d(:,:,imem,Nmon+it)   ! Luu vao Pr_All voi time tu Nmon+1 = it1 (Tinh tu   1/1983)
               Enddo
            Enddo
            Deallocate (xlon,xlat,xmem,xtime,X4d)
         Enddo
         Print*,"Reading data ...LT=",ilead
      Enddo
      Print*,"Nt=",Nt
      if (Allocated(xtime)) Deallocate (xtime) 
      Allocate (xtime(Ntime))
      !pause;print*,time
      xtime(:)=time(:)   ! restore xtime
      Deallocate (time)
      Print*,"Reading data ... Done"
   End Subroutine Read_all_data
   !------------------------------------------------------------------------------
   Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
   Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"

   Include "/home/tanpv/Mylibs/Sub_SPEI_One_Stn_2023.f90"
   Include "/home/tanpv/Mylibs/Sub_SPI_One_Stn_2023.f90"

   !Include "/home/tanpv/Mylibs/Sub_SPI.f90"  ! Chuong trinh tinh chi so han SPI
   !Include "/home/tanpv/Mylibs/Sub_SPEI.f90" ! Chuong trinh tinh chi so han SPEI
   Include "/home/tanpv/Mylibs/No_missing_Ar1.f90"

   !--------------------------------------------------------------------
   END ! Program