Program Dr_IDXs_SEACLID
  ! Calculating PDSI, SPI, SPEI 
  use netcdf

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

  Real, parameter    :: Miss=-99
  Integer, Parameter :: Nmon = 12
  Integer, Parameter :: Nreg = 7, Nmem = 7, Nlead = 6
  Integer, Parameter :: Ncase = 3       !  mod,cli,ann
  !Character (Len=3), Parameter :: Cases(Ncase)=(/"mod","cli","ann","mlr"/)
  Character (Len=3), Parameter :: Cases(Ncase)=(/"cli","ann","mlr"/)

  Integer, Allocatable :: Len_SPI(:)
  Integer, Parameter :: NDr_Lev = 3
  Character (Len=4), Parameter :: Dr_levnm(NDr_Lev)=(/"mild","mord","seve"/)
  Real, Parameter    :: Thes(NDr_Lev)=(/-0.5, -1.0, -1.5/)
  Character (Len=5)    :: IDX_nm             ! Name of Drought Index

  Real, Allocatable :: Regs(:,:), Lon(:), Lat(:), xlon(:), xlat(:) 
  Real, Allocatable :: Dr_IDX_All(:,:,:,:,:)                  ! Nlon,Nlat,Nspl,Nlead,NLen
  Real, Allocatable :: Dr_IDX_ens(:,:,:,:)                   ! Nlon,Nlat,NLen,Nlead       - Ensemble mean/median of Dr_IDX
  Real, Allocatable :: Dr_prob(:,:,:,:,:)                     ! Nlon,Nlat,Nlead,NDr_Lev,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(:), IDXname(:)
  Logical :: status, Empty

  Integer :: Yr_Ope, Mon_Ope, Nmodel, NYear, maxyrs, Nidx
  Integer :: NLen, Nst, Nlon, Nlat, IO, Ntime, Nt, it, Landpoint
  Integer :: imod, ist, iyr, imw, ivar, isys, system, ilon, ilat
  Integer :: Nt1, Nspl 

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

  Namelist /File_Dir_Names/ Inp_Dir, Mask_File, Out_Dir
  Namelist /Data_Info/      Yr_Ope, Mon_Ope, NLen, Nidx, Nmodel
  Namelist /IDX_name/       IDXname
  Namelist /Dr_Timescale/   Len_SPI
!------------------------------------------------------------------------------

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

  Read (1, File_Dir_Names)
  Read (1, Data_Info)      ! Yr_Ope, Mon_Ope, NLen, Nidx, Nmodel
  Allocate (Len_SPI(NLen)); Allocate (IDXname(Nidx))
  Read (1, IDX_name)       ! "spi", "spei",
  Read (1, Dr_Timescale)   ! 1,3,6,
  Close (1)

  Nspl = Nmodel*Nmem*Ncase*Nidx  ! = 2*7*3*2
  !Ntime = (Yr2-Yr1+1)*Nmon; maxyrs = Yr2-Yr1+1; NYear = (Yr_Ope - Yr1 + 1)
  Call Read_Regs         ! Mask file (Regs)
  print*,"Parameters: Nlat,Nlon,Nst,Nlen,Ntime=",Nlat,Nlon,Nst,Nlen,Ntime

  Call Read_all_data     ! Read DrIDXx for case = 56 members = 2 models x 7 mems x 4 cases ("mod","cli","ann","mlr")
  Call Calc_Dr_prob      ! Probability of occurrence of drought levels = Number of cases among 56 members divide to 56 
  Call Calc_Dr_IDX_Ens   ! Ensemble mean/median of drought indices
  Call Write_Dr_Prob
   CONTAINS
   !------------------------------------------------------------------------------
   Subroutine Calc_Dr_prob
     Integer :: ilon,ilat,ilead,ilen,ithre,ispl, ic, isum
     Real :: tmp, ntmp

     ! Dr_prob(:,:,:,:,:)                     ! Nlon,Nlat,Nlead,NDr_Lev,NLen 
     ! Dr_IDX_All(:,:,:,:,:)                  ! Nlon,Nlat,Nspl,Nlead,NLen
     ! Thes(NDr_Lev)=(/-0.5, -1.0, -1.5/)
     Allocate (Dr_prob(Nlon,Nlat,Nlead,NDr_Lev,NLen))
     Dr_prob=Miss
     Do ilen=1,NLen
        Do ithre=1,NDr_Lev
           Do ilead=1,Nlead
              Do ilon=1,Nlon
                 Do ilat=1,Nlat
                    if (Int(Regs(ilon,ilat))>0) then
                       isum=0; ic=0
                       Do ispl=1,Nspl
                          if (Dr_IDX_All(ilon,ilat,ispl,ilead,iLen)/=Miss) then
                             isum=isum+1
                             if (Dr_IDX_All(ilon,ilat,ispl,ilead,iLen)<=Thes(ithre)) ic=ic+1
                          endif
                       Enddo
                       if (isum>0) Dr_prob(ilon,ilat,ilead,ithre,ilen)=real(ic)/real(isum)*100.0
!print*,ic,isum,Dr_prob(ilon,ilat,ilead,ithre,ilen)
                    endif
                 Enddo
              Enddo
           Enddo
        Enddo
     Enddo
   End Subroutine Calc_Dr_prob
   !------------------------------------------------------------------------------
   Subroutine Calc_Dr_IDX_Ens
      Integer :: ilon,ilat,ilead,ilen,ithre,ispl, ic, isum
      Real :: xtmp(Nspl), mean, median
   ! Dr_IDX_ens((:,:,:,:)                   ! Nlon,Nlat,NLen,Nlead 

      Allocate (Dr_IDX_ens(Nlon,Nlat,NLen,Nlead)); Dr_IDX_ens=Miss
      Do ilead=1,Nlead
         Do ilen=1,NLen
            Do ilon=1,Nlon
               Do ilat=1,Nlat
                  if (Int(Regs(ilon,ilat))>0) then
                     xtmp(1:Nspl)=Dr_IDX_All(ilon,ilat,1:Nspl,ilead,iLen)
                     mean   = Ave (xtmp,Nspl,Miss) 
                     median = percentile (xtmp,Nspl,0.5,Miss)
                     !Dr_IDX_ens(ilon,ilat,ilen,ilead)=percentile (xtmp,Nspl,0.5,Miss)
                     if (mean/=Miss.AND.median/=Miss) then
                        Dr_IDX_ens(ilon,ilat,ilen,ilead) = median !(mean+median)/2.
                     else
                        Dr_IDX_ens(ilon,ilat,ilen,ilead) = Miss
                     endif
!print"(5F8.3)",mean,median,Dr_IDX_ens(ilon,ilat,ilen,ilead)
                  endif
               Enddo
            Enddo
         Enddo
      Enddo
   End Subroutine Calc_Dr_IDX_Ens
   !------------------------------------------------------------------------------
   Subroutine Write_Dr_Prob !(IDX_nm,ilen)
     !Character (Len=*) :: IDX_nm
     Character (Len=250) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm
     Integer :: i,ilen,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet
     Real :: tmpr, tmpt, tmpet, alpha, lat_reg
     Real :: xlead(Nlead)
     Real, Allocatable :: X4d(:,:,:,:),xlen(:)

     ! Write out for testing: SPEI_Out(ilon,ilat,imem,ilead,icas,imod,ilen)
     ! X4d(Nlon,Nlat,Nz,Nlead): Nz=Nmodel*Nmem*Ncase

     ! Dr_prob(:,:,:,:,:)                     ! Nlon,Nlat,Nlead,NDr_Lev,NLen

     write (Fname,"(A,I2.2,A)") trim(Out_Dir)//"Dr_Prob.nc"

     Allocate (X4d(Nlon,Nlat,Nlen,Nlead),xlen(Nlen))

     write(t_unit,"(A,I4.4,A,I2.2,A)") "months since ",Yr_Ope,"-",Mon_Ope,"-01 " 
     var_unit="%";
     lon_nm="lon"; lat_nm="lat"; z_nm="timescale"; t_nm="lead"
     Do ilead=1,Nlead; xlead(ilead)=ilead; Enddo
     Do ilen=1,Nlen; xlen(ilen)=Len_SPI(ilen); Enddo

     Do i=1,NDr_Lev
        Do ilen=1,Nlen
           Do ilead=1,Nlead
              X4d(:,:,ilen,ilead)=Dr_prob(:,:,ilead,i,ilen)
!do ilat=Nlat-2,Nlat-8,-1
!   print"(2I3,15F7.2)",ilen,ilead,X4d(14:20,ilat,ilen,ilead)
!enddo
           Enddo
        Enddo
!pause
        Var_ln="Probability of "//trim(Dr_levnm(i))//" drought occurrence. z=1..3: Timescale 1,3,6 months. t=1..6: Leadtimes"
        if (i==1) then
           Call wr_ncdf_4d (Fname,Lon,Lat,xlen,xlead, Nlon,Nlat,Nlen,Nlead,X4d,trim(Dr_levnm(i)), &
                              t_unit,Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss)
        else
           Call wr_ncdf_4d_apd (Fname,X4d,trim(Dr_levnm(i)),Nlon,Nlat,Nlen,Nlead, &
                              Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss)
        endif
     Enddo

     ! Ensemble mean/median of drought indices (Dr_IDX_ens):
     Var_ln="Ensemble mean/median of drought indices. z=1..3: Timescale 1,3,6 months. t=1..6: Leadtimes"
     Call wr_ncdf_4d_apd (Fname,Dr_IDX_ens,"dr_ens",Nlon,Nlat,Nlen,Nlead, &
                              Var_ln,"-",lon_nm, lat_nm, z_nm, t_nm,Miss)

     !Call wr_ncdf_4d (Fname,Lon,Lat,xcase,xlead, Nlon,Nlat,Nz,Nlead,X4d,"dridx", &
     !                         t_unit,Var_ln,"-",lon_nm, lat_nm, z_nm, t_nm,Miss)
     !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)
     !-----------------------------------------
     Deallocate (X4d,xlen)
   End Subroutine Write_Dr_Prob
   !------------------------------------------------------------------------------
   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,Lon,Lat,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 :: ispl,ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nle,icas, it1, it2
      Character (Len=255) :: Fname
      Real, Allocatable :: X4d(:,:,:,:), xlead(:), xlon(:),xlat(:),xmem(:) !,xtime(:)

      ! Nspl = Nmodel*Nmem*Ncase*Nidx
      ! Dr_IDX_All(:,:,:,:,:)                  ! Nlon,Nlat,Nspl,Nlead,NLen 
      Allocate (Dr_IDX_All(Nlon,Nlat,Nspl,Nlead,Nlen))
      Dr_IDX_All=Miss
      if (Allocated(X4d))  Deallocate (X4d);  if (Allocated(xmem))  Deallocate (xmem)
      if (Allocated(xlon)) Deallocate (xlon); if (Allocated(xlat))  Deallocate (xlat)
      
      Do ilen=1,Nlen
         ispl=0
         Do iidx=1,Nidx
            write (Fname,"(A,I2.2,A)") trim(Inp_Dir)//"Dr_IDX_"//trim(IDXname(iidx))//"_",Len_SPI(ilen),".nc"
            print*,"Reading ",trim(Fname)
            Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xlead,Nlo,Nla,Nme,Nle,X4d,"dridx")
!            Where (X4d>10.) X4d=Miss
!do ilon=1,Nlo
!   do ilat=1,Nla
!      do ie=1,Nme
!         do il=1,Nle
!            if (X4d(ilon,ilat,ie,il)>10) print*,ilon,ilat,ie,il,X4d(ilon,ilat,ie,il);
!         enddo
!      enddo
!   enddo
!enddo
!pause
            Do imem=1,Nme
               ispl=ispl+1
               Do ilead=1,Nle
                  Dr_IDX_All(:,:,ispl,ilead,ilen)=X4d(:,:,imem,ilead) 
               Enddo
               !print*,"ilen,ispl=",ilen,ispl
            Enddo
            Deallocate (xlon,xlat,xmem,xlead,X4d)
         Enddo
      Enddo
      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/Ave.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