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    :: beta_spi(:), gamm(:), pzero(:), SPI(:,:,:), SPE(:,:,:)
  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,Nmem,Ntime,Ncase,Nlen

  ! 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

  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

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

  Do imod=1,1 !Nmodel ==> Cho OBS chi can 1 truong hop
     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

     ! For calculating SPI, SPEI:
     ! Do voi OBS tat ca cac model, cac mem va cac lead deu nhu nhau ==> Chi lay 1 tr/h
      Ntime1=Ntime           !!!!!!!!!!
      Do ilead=1,1 !Nlead
         Nt1 = Nt !+ ilead    ! Nam cuoi cung cac leadtime co SL keo dai them. Doi voi OBS thi khong can
         SPEI_All=Miss
         Do imem=1,1 !Nmem
            Do icas=1,Ncase   ! Ncase=1=OBS
               ist=0; DatBL=Miss; DatTMP=Miss; Ntime=Ntime1
               Do ilon=1,Nlon
                  Do ilat=1,Nlat
                     if (Regs(ilon,ilat)>0) then
                        ist=ist+1;
                        DatBL (1:Ntime,ist) = Pr_All(ilon,ilat,imem,ilead,1:Ntime,icas) !Rain_Inp(ilon,ilat,1:Ntime)
                        DatTMP(1:Ntime,ist) = T2_All(ilon,ilat,imem,ilead,1:Ntime,icas) !Temp_Inp(ilon,ilat,1:Ntime)
                        Lat(ist)=xlat(ilat)
                        !X1d(1:Ntime)=DatBL (1:Ntime,ist); Nt1=Length_series (X1d,Ntime,Miss); print"(A,8I6)","ilead,imem,icas: ",ilead,imem,icas,Nt1 
                        !X1d(1:Ntime)=DatTMP(1:Ntime,ist); Nt1=Length_series (X1d,Ntime,Miss); print"(A,8I6)","ilead,imem,icas: ",ilead,imem,icas,Nt1
                        !X1d(1:Ntime)=DatTMP(1:Ntime,ist); Nt1=Length_series (X1d,Ntime,Miss)
                        !if (Nt1>Ntime1/3) then; Ntime=Nt1; Nt=Nt1; endif
                     endif
                  Enddo
               Enddo
!pause
               Ntime=Nt1
               if (Allocated(SPI)) Deallocate (SPI); Allocate (SPI(Nst,Ntime,NLen))
               SPI(:,:,:) = amssng
               Call Cal_SPEI
               !if (icas==1.AND.imem==3) then; print*,SPI,"LT=",ilead; pause; endif
               ist=0
               Do ilon=1,Nlon
                  Do ilat=1,Nlat
                     if (Regs(ilon,ilat)>0) then
                        ist=ist+1
                        Do ilen=1,Nlen
                           SPEI_All(ilon,ilat,imem,1:Ntime,icas,ilen)=SPI(ist,1:Ntime,ilen)
                        Enddo
                        !X1d(1:Ntime)=SPI(ist,1:Ntime,1); Nt1=Length_series (X1d,Ntime,Miss); print"(A,8I6)","ilead,imem,icas: ",ilead,imem,icas,Nt1
                     endif
                  Enddo
               Enddo
!pause
            Enddo
         Enddo
         !iyr=Nt1/Nmon; imon=Nt1-iyr*Nmon; iyr=Yr1+iyr-1;
         !print*,"Years: ",Yr1,"-",iyr,"/",imon; pause
         Ntime=Ntime1
         !write (Fname,"(A,I1,A)") trim(Out_Dir)//"SPEI_"//trim(Model_nm(imod))//"_lead_",ilead,".nc"
         write (Fname,"(A,I4.4,A,I4.4,A)") trim(Out_Dir)//"SPEI_obs",Yr1,"-",Yr2,".nc"
         Call Write_Dr_IDXs_to_file (Fname,ilead,Model_nm(imod)) 
      Enddo
   Enddo
   CONTAINS
   !------------------------------------------------------------------------------
   Subroutine Reinit
      if (Allocated(T2_All)) Deallocate (T2_All)
      if (Allocated(Pr_All)) Deallocate (Pr_All)

      if (Allocated (Rain))     Deallocate (Rain)
      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 (beta_spi(Nmonth), gamm(Nmonth), pzero(Nmonth))
      !Allocate (SPI(Nst,Ntime,NLen), SPE(Nst,Ntime,NLen))
      Allocate (SPEI_All(Nlon,Nlat,Nmem,Ntime,Ncase,Nlen))
      Allocate (Lat(Nst),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, thres, xxx(Ntime)
     Real, Allocatable :: X41(:,:,:,:)

     t_unit="months since 1982-01-01 "; var_unit="-";
     lon_nm="lon"; lat_nm="lat"; z_nm="timescale"; t_nm="time"
     Do imem=1,Nlen; xlen(imem)=Len_SPI(imem);  Enddo
     Allocate (X4d(Nlon,Nlat,Nlen,Ntime))
     !Do icas=1,Ncase
        Do ilen=1,Nlen 
           !Do imem=1,1 !Nmem 
              Do ilon=1,Nlon
                 Do ilat=1,Nlat
                    Do it=1,Ntime
                       X4d(ilon,ilat,ilen,it)=SPEI_All(ilon,ilat,1,it,1,ilen)
                    Enddo
                 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" 
           !write(Var_nm,"(A,I2.2)") trim(Cases(icas))//"_",Len_SPI(ilen)
           !if (icas==1.AND.ilen==1) then
        Call wr_ncdf_4d (Fname, xlon,xlat,xlen,xtime, Nlon,Nlat,Nlen,Ntime,X4d,"idx", &
                              t_unit,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));
        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)
                    if (tmp/=Miss.AND.tmp<=thres) X41(ilon,ilat,ilen,it)=1
                    if (tmp/=Miss.AND.tmp>thres)  X41(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,"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)
                    if (tmp/=Miss.AND.tmp<=thres) X41(ilon,ilat,ilen,it)=1
                    if (tmp/=Miss.AND.tmp>thres)  X41(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,"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)
                    if (tmp/=Miss.AND.tmp<=thres) X41(ilon,ilat,ilen,it)=1
                    if (tmp/=Miss.AND.tmp>thres)  X41(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,"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)
                    if (tmp/=Miss.AND.tmp<=thres) X41(ilon,ilat,ilen,it)=1
                    if (tmp/=Miss.AND.tmp>thres)  X41(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,"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)
   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=1,Ninp !,1,-1
        if (X(i)/=Miss) then
           j=j+1
        else
           exit
        endif
     Enddo
     !print*,Ninp,j
     Length_series=j !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_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