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 = 4       !  obs,mod,cli,ann
  Integer, Parameter :: Ncase_1 = 2     ! obs, ens
  Character (Len=3), Parameter :: Case_1(Ncase_1)=(/"obs","ens"/)
  Character (Len=3), Parameter :: Cases(Ncase)=(/"obs","mod","cli","ann"/)
  !Character (Len=3), Parameter :: Cases(Ncase)=(/"obs","cli","ann"/)

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

  ! Distribution function's parameters:
  Real, Allocatable    :: beta_spi(:), gamm(:), pzero(:), SPI(:,:,:), SPE(:,:,:)
  Real                 :: beta_spei(0:2), logLogisticParams(0:Nmon,0:2), DisPars(0:Nmon,0:2)
  Logical              :: use_baseline
  Character (Len=5)    :: IDX_nm             ! Name of Drought Index

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

!------------------------------------------------------------------------------
  ! Ensemble:
  Real, Allocatable :: Pr_ens  (:,:,:,:,:), T2_ens(:,:,:,:,:)       ! Nlon,Nlat,Nlead,Ntime,Ncase_1
  Real, Allocatable :: Pet_ens (:,:,:,:,:)
  Real, Allocatable :: SPEI_ens(:,:,:,:,:)                          ! Nlon,Nlat,Ntime,Ncase_1,Nlen

  Real, Allocatable :: Pr_Reg (:,:,:,:), T2_Reg(:,:,:,:) ! (Nreg,Nlead,Ncase_1,Ntime) - Obs va Ens
  Real, Allocatable :: Pet_Reg(:,:,:,:)
  Real, Allocatable :: SPI_Reg(:,:,:,:,:)                ! (Nreg,Nlead,Ncase_1,Ntime,Nlen)

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

  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)

  Nyear = (Yr2 - Yr1 + 1); Ntime = Nyear*NMonth; maxyrs = NYear
  Nt = Ntime - Nmon   ! 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

  Call Read_all_data
  Call Do_Ensemble
  Call Reinit
  Call Calc_IDXs ("spei") 
  Call Areal_Mean           ! Regional means of Pr, T2, Pet
  Call Calc_Regional_Dr_IDXs ("spei")
  Call Write_Pr_Temp_DrIDXs ("spei")
   CONTAINS
   !------------------------------------------------------------------------------
   Subroutine Calc_IDXs (IDX_nm)   ! IDX_nm="spi", "spei"
      Character (Len=*) :: IDX_nm
      Integer :: ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet
      Real :: tmpr, tmpt, tmpet, alpha, latx
      Real :: xspi(Ntime)

     ! For calculating SPI, SPEI:
      SPEI_ens=Miss
      Do ilead=1,Nlead
         Do ilen=1,Nlen
            Do icas=1,Ncase_1
               Do ilon=1,Nlon
                  Do ilat=1,Nlat
                     if (Regs(ilon,ilat)>0) then
                        use_baseline=.false.; xspi=Miss; latx=xlat(ilat)
                        Rain(1:Ntime)=Pr_ens(ilon,ilat,ilead,1:Ntime,icas)
                        T2m (1:Ntime)=T2_ens(ilon,ilat,ilead,1:Ntime,icas)
                        Nt1=Min(Length_series (Rain,Ntime,Miss), Length_series (T2m,Ntime,Miss))
                        if (IDX_nm=="spi") then
                           Call Cal_SPI_One_Stn (Rain,Nt1,Len_SPI(ilen),alpha,beta_spi,gamm,pzero,xspi,use_baseline,Miss)
                           SPEI_ens(ilon,ilat,1:Nt1,icas,ilen)=xspi(1:Nt1)
                        else if (IDX_nm=="spei") then
                           Call Calc_SPEI_One_Stn (T2m,Rain,latx,Nt1,Nmon,Len_SPI(ilen),logLogisticParams,use_baseline,xspi,Miss)
                           SPEI_ens(ilon,ilat,1:Nt1,icas,ilen)=xspi(1:Nt1)
                        endif
                     endif
                  Enddo
               Enddo
            Enddo
         Enddo
         write (Fname,"(A,I1,A)") trim(Out_Dir)//trim(IDX_nm)//"_ens_lead_",ilead,".nc"
         Call Write_Dr_IDXs_to_file (Fname,ilead) !,Model_nm(imod)) 
      Enddo
   End Subroutine Calc_IDXs
   !------------------------------------------------------------------------------
   Subroutine Calc_Regional_Dr_IDXs (IDX_nm)   ! IDX_nm="spi", "spei"
     Character (Len=*) :: IDX_nm
     Integer :: ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet
     Real :: tmpr, tmpt, tmpet, alpha, lat_reg
     Real :: xspi(Ntime)
     Logical :: use_baseline    ! =.false. if distribution's parameters (alpha, beta,...) need to be calculated
                                ! =.true.  if not

     SPI_Reg = Miss; use_baseline=.false.
     Do ilen=1,Nlen
        Do ireg=1,Nreg
           Do ilead=1,Nlead
              Do icas=1,Ncase_1
                    Rain(1:Ntime)=Pr_Reg (ireg,ilead,icas,1:Ntime)
                    T2m (1:Ntime)=T2_Reg (ireg,ilead,icas,1:Ntime)
                    PET (1:Ntime)=Pet_Reg(ireg,ilead,icas,1:Ntime)
                    Nt = Min(Length_series (Rain,Ntime,Miss), Length_series (T2m,Ntime,Miss))
                    if (IDX_nm=="spi") then
                       Call Cal_SPI_One_Stn (Rain,Nt,Len_SPI(ilen),alpha,beta_spi,gamm,pzero,xspi,use_baseline,Miss)
                       SPI_Reg(ireg,ilead,icas,1:Nt,ilen)=xspi(1:Nt)
                    else if (IDX_nm=="spei") then
                       print"(24F6.1)",Rain(1:Ntime)
                       print*,"icas, Nt,use_baseline=",icas,Cases(icas),Nt,use_baseline
                       Call Calc_SPEI_One_Stn_using_PET (Rain, PET, Nt, Nmon, Len_SPI(ilen), logLogisticParams, use_baseline, xspi, Miss)
                       DisPars = logLogisticParams
                       SPI_Reg(ireg,ilead,icas,1:Nt,ilen)=xspi(1:Nt)
                    endif
              Enddo
           Enddo
        Enddo
     Enddo
   End SUbroutine Calc_Regional_Dr_IDXs
   !------------------------------------------------------------------------------
   Subroutine Write_Pr_Temp_DrIDXs (IDX_nm)
     Character (Len=*) :: IDX_nm
     Character (Len=250) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm
     Integer :: ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet
     Real :: tmpr, tmpt, tmpet, alpha, lat_reg
     Real :: xreg(Nreg),xlead(Nlead),xcase(Ncase_1),xspi(Ntime)
     Real, Allocatable :: X4d(:,:,:,:)

     ! Write out for testing:
     !Fname="Areal_mean_Pr_T2_"//trim(IDX_nm)//"_"//trim(Model_nm)//".nc"
     write (Fname,"(A,I1,A)") trim(Out_Dir)//"Areal_mean_"//trim(IDX_nm)//"_ens.nc"

     t_unit="months since 1983-01-01 "; var_unit="-";
     lon_nm="reg"; lat_nm="lead"; z_nm="case"; t_nm="time"
     Do ireg=1,Nreg; xreg(ireg)=ireg; Enddo
     Do ilead=1,Nlead; xlead(ilead)=ilead; Enddo
     Do icas=1,Ncase_1; xcase(icas)=icas; Enddo
     Var_ln="Areal mean of Precipitation: Z=1..Ncase_1: obs,ens"
     Call wr_ncdf_4d (Fname, xreg,xlead,xcase,xtime, Nreg,Nlead,Ncase_1,Ntime,Pr_Reg,"pr", &
                              t_unit,Var_ln,"mm/month",lon_nm, lat_nm, z_nm, t_nm,Miss)
     Var_ln="Areal mean of T2m: Z=1..Ncase_1: obs,ens"
     Call wr_ncdf_4d_apd (Fname,T2_Reg,"t2",Nreg,Nlead,Ncase_1,Ntime, &
                              Var_ln,"C",lon_nm, lat_nm, z_nm, t_nm,Miss)
     Var_ln="Areal mean of PET: Z=1..Ncase_1: obs,ens"
     Call wr_ncdf_4d_apd (Fname,Pet_Reg,"pet",Nreg,Nlead,Ncase_1,Ntime, &
                              Var_ln,"mm/month",lon_nm, lat_nm, z_nm, t_nm,Miss)
     !-----------------------------------------
     Allocate (X4d(Nreg,Nlead,Ncase_1,Ntime))
     Do ilen=1,Nlen
        X4d(:,:,:,:)=SPI_Reg(:,:,:,:,ilen)
        write(Var_nm,"(A,I1)") "spi",Len_SPI(ilen)
        write(Var_ln,"(A,I1)") "SPI_",Len_SPI(ilen)
        Call wr_ncdf_4d_apd (Fname,X4d,trim(Var_nm),Nreg,Nlead,Ncase_1,Ntime, &
                              Var_ln,"-",lon_nm, lat_nm, z_nm, t_nm,Miss)
     Enddo
   End Subroutine Write_Pr_Temp_DrIDXs
   !------------------------------------------------------------------------------
   Subroutine Areal_Mean !(Model_nm)
     Integer :: ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet
     Integer :: Ngrd, ic1, ic2, ic3
     Real :: tmpr, tmpt, tmpet, alpha, lat_reg
     Real, Allocatable, Dimension (:) :: xtmp1, xtmp2, xtmp3

     Pr_Reg=Miss; T2_Reg=Miss; Pet_Reg=Miss
     Do ilead=1,Nlead
        Do icas=1,Ncase_1
           Do it=1,Ntime
              Do ireg=1,Nreg
                 Ngrd = Get_Reg_Points (ireg)
                 Allocate (xtmp1(Ngrd), xtmp2(Ngrd), xtmp3(Ngrd))
                 ic1=0; ic2=0; ic3=0
                 Do ilon=1,Nlon
                    Do ilat=1,Nlat
                       if (Int(Regs(ilon,ilat))==ireg) then
                          !if (Pr_All(ilon,ilat,imem,ilead,it,icas)/=Miss) then
                          if (Pr_ens(ilon,ilat,ilead,it,icas)/=Miss) then
                             ic1=ic1+1; xtmp1(ic1)=Pr_ens(ilon,ilat,ilead,it,icas)
                          endif
                          if (T2_ens(ilon,ilat,ilead,it,icas)/=Miss) then
                             ic2=ic2+1; xtmp2(ic2)=T2_ens(ilon,ilat,ilead,it,icas)
                          endif
                          if (Pet_ens(ilon,ilat,ilead,it,icas)/=Miss) then
                             ic3=ic3+1; xtmp3(ic2)=Pet_ens(ilon,ilat,ilead,it,icas)
                          endif
                       endif
                    Enddo
                 Enddo
                 Pr_Reg (ireg,ilead,icas,it)=Ave(xtmp1,ic1,Miss) !percentile (xtmp1,ic1,0.5,Miss)
                 T2_Reg (ireg,ilead,icas,it)=Ave(xtmp2,ic2,Miss) !percentile (xtmp2,ic2,0.5,Miss)
                 Pet_Reg(ireg,ilead,icas,it)=Ave(xtmp3,ic3,Miss) !percentile (xtmp3,ic3,0.5,Miss)
                 Deallocate (xtmp1,xtmp2,xtmp3)
              Enddo
           Enddo
        Enddo
     Enddo
   End Subroutine Areal_Mean
   !------------------------------------------------------------------------------
   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 (PET))      Deallocate (PET)

      if (Allocated (xspi))     Deallocate (xspi)

      if (Allocated (Pr_Reg))   Deallocate (Pr_Reg)
      if (Allocated (T2_Reg))   Deallocate (T2_Reg)
      if (Allocated (Pet_Reg))  Deallocate (Pet_Reg)
      if (Allocated (beta_spi)) Deallocate (beta_spi)
      if (Allocated (gamm))     Deallocate (gamm)
      if (Allocated (pzero))    Deallocate (pzero)
      if (Allocated (SPEI_ens)) Deallocate (SPEI_ens)
      if (Allocated (Lat))      Deallocate (Lat)
      if (Allocated (X1d))      Deallocate (X1d)

      Allocate (Rain(Ntime),T2m(Ntime),PET(Ntime),xspi(Ntime)) !DatBL(Ntime,Nst),DatTMP(Ntime,Nst))
      Allocate (beta_spi(Nmonth), gamm(Nmonth), pzero(Nmonth))
      Allocate (SPEI_ens(Nlon,Nlat,Ntime,Ncase_1,Nlen))                          ! Nlon,Nlat,Ntime,Ncase_1,Nlen

      Allocate (Pr_Reg(Nreg,Nlead,Ncase_1,Ntime),T2_Reg(Nreg,Nlead,Ncase_1,Ntime))
      Allocate (Pet_Reg(Nreg,Nlead,Ncase_1,Ntime))
      Allocate (SPI_Reg(Nreg,Nlead,Ncase_1,Ntime,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 :: zlev(Nlen)
     Real, Allocatable :: X3d(:,:,:)

     t_unit="months since 1983-01-01 "; var_unit="-";
     lon_nm="lon"; lat_nm="lat"; z_nm="timescale"; t_nm="time"
     Do ilen=1,Nlen; zlev(ilen)=Len_SPI(ilen);  Enddo
     Allocate (X4d(Nlon,Nlat,Nlen,Ntime))
     ! SPEI_ens(Nlon,Nlat,Ntime,Ncase_1,Nlen)
     ! Pr_ens  (Nlon,Nlat,Nlead,Ntime,Ncase_1)
     Do icas=1,Ncase_1
        Do ilen=1,Nlen 
           Do ilon=1,Nlon
              Do ilat=1,Nlat
                 Do it=1,Ntime
                    X4d(ilon,ilat,ilen,it)=SPEI_ens(ilon,ilat,it,icas,ilen)
                 Enddo
              Enddo
           Enddo
        Enddo
        
        write(Var_ln,"(A,I2.2,A,I1)") "SPEI for "//trim(Case_1(icas))//" with timescale in z=1,2,3 corresponding to 1,3,6 month(s) for Leadtime=",ilead
        write(Var_nm,"(A,I2.2)") trim(Case_1(icas)) !//"_",Len_SPI(ilen)
        if (icas==1) then
           Call wr_ncdf_4d (Fname, xlon,xlat,zlev,xtime, Nlon,Nlat,Nlen,Ntime,X4d,Var_nm, &
                              t_unit,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,Nlen,Ntime, &
                              Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss)
        endif
     Enddo
     Deallocate (X4d)
     ! Write out Pr and T2m:
     Allocate (X3d(Nlon,Nlat,Ntime))
     X3d(:,:,1:Ntime)=Pr_ens (:,:,ilead,1:Ntime,1)    ! Pr_obs
     Call wr_ncdf_3d_apd (Fname,X3d,"pr_obs",Nlon,Nlat,Ntime,"Pr obs","mm/month",lon_nm,lat_nm,t_nm,Miss)
     X3d(:,:,1:Ntime)=Pr_ens (:,:,ilead,1:Ntime,2)    ! Pr_model_ens
     Call wr_ncdf_3d_apd (Fname,X3d,"pr_ens",Nlon,Nlat,Ntime,"Pr model ens","mm/month",lon_nm,lat_nm,t_nm,Miss)

     X3d(:,:,1:Ntime)=T2_ens (:,:,ilead,1:Ntime,1)    ! T2_obs
     Call wr_ncdf_3d_apd (Fname,X3d,"t2_obs",Nlon,Nlat,Ntime,"T2 obs","C",lon_nm,lat_nm,t_nm,Miss)
     X3d(:,:,1:Ntime)=T2_ens (:,:,ilead,1:Ntime,2)    ! Pr_model_ens
     Call wr_ncdf_3d_apd (Fname,X3d,"t2_ens",Nlon,Nlat,Ntime,"T2 model ens","C",lon_nm,lat_nm,t_nm,Miss)
     Deallocate (X3d)
   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
   !--------------------------------------------------------------------------------
   Integer Function Get_Reg_Points (ireg)
      Integer :: ireg,ilon,ilat,Npoint

      Npoint=0
      Do ilon=1,Nlon
         Do ilat=1,Nlat
            if (Regs(ilon,ilat)==ireg) Npoint=Npoint+1
         Enddo
      Enddo
      Get_Reg_Points = Npoint
   End Function Get_Reg_Points
   !------------------------------------------------------------------------------
   Subroutine Read_all_data ! (Model_nm)
      !Character (Len=*) :: Model_nm
      Integer :: imod,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 imod=1,Nmodel
         Do ilead=1,Nlead
            write (Fname,"(A,I2.2,A)") trim(Inp_Dir)//trim(TVar)//"_"//trim(Model_nm(imod))//"_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 (imod==1.AND.icas==1.AND.ilead==1) then 
                  Allocate (T2_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase,Nmodel),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,imod)=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(imod))//"_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 (imod==1.AND.icas==1.AND.ilead==1) then; Allocate (Pr_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase,Nmodel)); Pr_All=Miss; endif
               Do imem=1,Nmem
                  Do it=1,Ntime
                     Pr_All(:,:,imem,ilead,it,icas,imod)=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
      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
   !------------------------------------------------------------------------------
   Subroutine Do_Ensemble
      Integer :: ilon,ilat,it,ilead,imem,icas, i,j,k, N
      Real, Allocatable :: X1d(:),Y1d(:)
      Real :: temp,rain, latx, xtmp(0:Ntime-1),xpet(0:Ntime-1)

      N = Nmodel*Nmem*(Ncase-1)
      Allocate (X1d(N),Y1d(N))
      Allocate (T2_ens(Nlon,Nlat,Nlead,Ntime,Ncase_1),Pr_ens(Nlon,Nlat,Nlead,Ntime,Ncase_1))
      Allocate (Pet_ens(Nlon,Nlat,Nlead,Ntime,Ncase_1))
      T2_ens=Miss; Pr_ens=Miss; Pet_ens=Miss
      Do ilead=1,Nlead           ! Obs data
         Do it=1,Ntime
            T2_ens(:,:,ilead,it,1)=T2_All(:,:,1,ilead,it,1,1)
            Pr_ens(:,:,ilead,it,1)=Pr_All(:,:,1,ilead,it,1,1)
         Enddo
      Enddo
      Do ilon=1,Nlon             ! Model data, including original and corrected
         Do ilat=1,Nlat
         if (Int(Regs(ilon,ilat))>0) then
            latx=xlat(ilat)
            Do ilead=1,Nlead     
               Do it=1,Ntime
                  i=0
                  Do imod=1,Nmodel
                     Do imem=1,Nmem
                        Do icas=3,Ncase   ! 2:"mod",3:"cli",4:"ann"
                           i=i+1
                           X1d(i)=T2_All(ilon,ilat,imem,ilead,it,icas,imod)
                           Y1d(i)=Pr_All(ilon,ilat,imem,ilead,it,icas,imod)
                        Enddo
                     Enddo
                  Enddo
                  temp=percentile (X1d,i,0.5,Miss); rain=percentile (Y1d,i,0.5,Miss)     ! Median instead of Mean
                  T2_ens(ilon,ilat,ilead,it,2)=temp
                  Pr_ens(ilon,ilat,ilead,it,2)=rain
               Enddo
               !
               xtmp(0:Ntime-1)=T2_ens(ilon,ilat,ilead,1:Ntime,1)    ! obs
               Call PET_thornthwaite (xtmp, Ntime, latx, xpet, Miss)
               Pet_ens(ilon,ilat,ilead,1:Ntime,1)=xpet(0:Ntime-1)
               xtmp(0:Ntime-1)=T2_ens(ilon,ilat,ilead,1:Ntime,2)    ! ens
               Call PET_thornthwaite (xtmp, Ntime, latx, xpet, Miss)
               Pet_ens(ilon,ilat,ilead,1:Ntime,2)=xpet(0:Ntime-1)
            Enddo
         endif
         Enddo
      Enddo
   End Subroutine Do_Ensemble
   !------------------------------------------------------------------------------
   Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
   Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
   Include "/home/tanpv/Mylibs/Ave.f90"

   Include "/home/tanpv/Mylibs/Sub_SPI_One_Stn_2023.f90"
   Include "/home/tanpv/Mylibs/Sub_SPEI_One_Stn_2023.f90"
   Include "/home/tanpv/Mylibs/PET_thornthwaite.f90"
   Include "/home/tanpv/Mylibs/Sub_SPEI_One_Stn_PET_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/Sub_Percentile.f90"
   !--------------------------------------------------------------------
   END ! Program