Program Dr_IDXs_SEACLID
  ! Calculating PDSI, SPI, SPEI 
  use netcdf

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

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

  Real, Allocatable :: Regs(:,:), Lon(:), Lat(:), xlon(:), xlat(:) 

  Real, Allocatable :: Pr_Ope(:,:,:,:,:,:), T2_Ope(:,:,:,:,:,:)     ! Nlon,Nlat,Nmem,Nlead,Ncase,Nmodel
  Real, Allocatable :: Pr_Ope_ens(:,:,:), T2_Ope_ens(:,:,:)       ! Nlon,Nlat,Nlead
  Real, Allocatable :: Pr_Ope_fcs(:,:,:), T2_Ope_fcs(:,:,:)       ! Nlon,Nlat,Nlead = *Cli_Obs + *Ano
  Real, Allocatable :: Pr_Ano(:,:,:), T2_Ano(:,:,:)               ! Nlon,Nlat,Nlead
  Real, Allocatable :: Pr_Cli(:,:,:,:), T2_Cli(:,:,:,:)           ! Nlon,Nlat,Nlead,Nmon
  Real, Allocatable :: Pr_Cli_Obs(:,:,:), T2_Cli_Obs(:,:,:)       ! Nlon,Nlat,Nmon

  Real, Allocatable :: xmem(:), xlead(:), xtime(:)

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

  Integer :: Yr1, Yr2, YOpe, MOpe, 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/ Cli_Dir, Mask_File, Dir_Ope, Dir_Out, TVar, RVar
  Namelist /Data_Info/      Nmodel
  Namelist /Model_name/     Model_nm
  Namelist /Time_Infor/     YOpe, MOpe
!------------------------------------------------------------------------------

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

  Read (1, File_Dir_Names)
  Read (1, Data_Info)         
  Allocate (Model_nm(Nmodel))
  Read (1, Model_name)
  Read (1, Time_Infor)
  Close (1)

  Call Read_Regs 
  print*,"Parameters: Nlat,Nlon=",Nlat,Nlon
  
  Call Read_Ope_Data
  Call Read_Clim_Data
  Call Do_ens_Ano_Ope
  Call Write_Forecast_Mean_Ano_Pr_T2 
   CONTAINS
   !------------------------------------------------------------------------------
   Subroutine Do_ens_Ano_Ope
      Integer :: imod,imem,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,icas, it1, it2
      Integer :: Nspl,islp
      Real, Allocatable :: xtmp(:), xpr(:)

      ! Pr_Cli(:,:,:,:), T2_Cli(:,:,:,:)           ! Nlon,Nlat,Nlead,Nmon
      ! Pr_Ope_fcs(:,:,:), T2_Ope_fcs(:,:,:)       ! Nlon,Nlat,Nlead
      ! Pr_Cli_Obs(:,:,:), T2_Cli_Obs(:,:,:)       ! Nlon,Nlat,Nmon

      Nspl=Nmodel*Nmem*Ncase
      Allocate (xtmp(Nspl),xpr(Nspl))
      Allocate (Pr_Ope_ens(Nlon,Nlat,Nlead),T2_Ope_ens(Nlon,Nlat,Nlead))
      Allocate (Pr_Ope_fcs(Nlon,Nlat,Nlead),T2_Ope_fcs(Nlon,Nlat,Nlead))
      Pr_Ope_ens=Miss; T2_Ope_ens=Miss; Pr_Ope_fcs=Miss; T2_Ope_fcs=Miss
      Allocate (Pr_Ano(Nlon,Nlat,Nlead),T2_Ano(Nlon,Nlat,Nlead))
      Pr_Ano=Miss; T2_Ano=Miss

      Do ilead=1,Nlead
         imon=MOpe+ilead; if (imon>Nmon) imon=imon-Nmon    ! Month of the year
         Do ilon=1,Nlon
            Do ilat=1,Nlat
               if (Int(Regs(ilon,ilat))>0) then
                  islp=0
                  Do imod=1,Nmodel
                     Do imem=1,Nmem
                        Do icas=1,Ncase
                           islp=islp+1
                           xpr (islp)=Pr_Ope(ilon,ilat,imem,ilead,icas,imod)
                           xtmp(islp)=T2_Ope(ilon,ilat,imem,ilead,icas,imod)
                        Enddo
                     Enddo
                  Enddo
                  ! Fcs ensemble mean:
                  Pr_Ope_ens(ilon,ilat,ilead)=Ave(xpr, islp,Miss)
                  T2_Ope_ens(ilon,ilat,ilead)=Ave(xtmp,islp,Miss)
                  ! Fcs Anomaly:
                  Pr_Ano(ilon,ilat,ilead)=Pr_Ope_ens(ilon,ilat,ilead)-Pr_Cli(ilon,ilat,ilead,imon) ! Ano = Ens_fcs - Mod_Cli
                  T2_Ano(ilon,ilat,ilead)=T2_Ope_ens(ilon,ilat,ilead)-T2_Cli(ilon,ilat,ilead,imon)
                  ! Operational fcs:
                  Pr_Ope_fcs(ilon,ilat,ilead)=Pr_Cli_Obs(ilon,ilat,imon)+Pr_Ano(ilon,ilat,ilead)   ! Fcs = Obs_Cli + Ano
                  T2_Ope_fcs(ilon,ilat,ilead)=T2_Cli_Obs(ilon,ilat,imon)+T2_Ano(ilon,ilat,ilead)
                  ! Relative anomaly for rainfall:
                  Pr_Ano(ilon,ilat,ilead)=Pr_Ano(ilon,ilat,ilead)/Pr_Cli_Obs(ilon,ilat,imon)       ! Relative Ano = Ano/Obs_Cli
               endif
            Enddo
         Enddo
      Enddo
   End Subroutine Do_ens_Ano_Ope
   !------------------------------------------------------------------------------
   Subroutine Write_Forecast_Mean_Ano_Pr_T2
      Character (Len=250) :: var_nm, t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm
      Real, Allocatable, Dimension (:) :: xtime
      INteger :: i, Nvar1

!      Allocate (Pr_Ope_ens(Nlon,Nlat,Nlead),T2_Ope_ens(Nlon,Nlat,Nlead))
!      Pr_Ope_ens=Miss; T2_Ope_ens=Miss
!      Allocate (Pr_Ano(Nlon,Nlat,Nlead),T2_Ano(Nlon,Nlat,Nlead))

      Allocate (xtime(Nlead))
      Do i=1,Nlead; xtime(i)=i; Enddo
      Write(t_unit,"(A,I4.4,A,I2.2,A)") "months since ",YOpe,"-",MOpe+1,"-01 "; 
      lon_nm="lon"; lat_nm="lat"; t_nm="lead";

      Fname=trim(Dir_Out)//"Forecast_Ope_Pr_T2_and_Anomaly.nc"
      var_ln="Model rainfall ensemble forecast"
      var_unit="mm/month"
      Call wr_ncdf_3d (Fname,Lon,Lat,xtime,Nlon,Nlat,Nlead,Pr_Ope_ens,"pr",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      var_ln="Model 2-m temperature ensemble forecast"
      var_unit="C degree"
      Call wr_ncdf_3d_apd (Fname,T2_Ope_ens,"t2",Nlon,Nlat,Nlead, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      var_ln="Relative Anomaly of the Model rainfall forecast"
      var_unit="-"
      Call wr_ncdf_3d_apd (Fname,Pr_Ano,"pr_ano",Nlon,Nlat,Nlead, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      var_ln="Anomaly of the Model 2-m temperature forecast"
      var_unit="C degree"
      Call wr_ncdf_3d_apd (Fname,T2_Ano,"t2_ano",Nlon,Nlat,Nlead, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)

      ! Pr_Ope_fcs(:,:,:), T2_Ope_fcs(:,:,:)       ! Nlon,Nlat,Nlead
      var_ln="Model rainfall forecast obtained by adding anomaly forecast to observation climatology"
      var_unit="mm/month"
      Call wr_ncdf_3d_apd (Fname,Pr_Ope_fcs,"pr_fcs",Nlon,Nlat,Nlead, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      var_ln="Model 2-m temperature forecast obtained by adding anomaly forecast to observation climatology"
      var_unit="C degree"
      Call wr_ncdf_3d_apd (Fname,T2_Ope_fcs,"t2_fcs",Nlon,Nlat,Nlead, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      Deallocate (xtime)
   End Subroutine Write_Forecast_Mean_Ano_Pr_T2
   !--------------------------------------------------------------------------------
   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
   !------------------------------------------------------------------------------
   Subroutine Read_Clim_Data
      Integer :: Nlo,Nla,Nle,Nt
      Real, Allocatable :: xlon(:),xlat(:),xled(:),xtime(:)

      ! Pr_Cli_Obs(:,:,:), T2_Cli_Obs(:,:,:)       ! Nlon,Nlat,Nmon

      Fname=trim(Cli_Dir)//"Monthly_Clim.nc"
      Call rd_ncdf_4d (Fname,xlon,xlat,xled,xtime,Nlo,Nla,Nle,Nt,Pr_Cli,"pr_m")
      Deallocate (xlon,xlat,xled,xtime)
      Call rd_ncdf_4d (Fname,xlon,xlat,xled,xtime,Nlo,Nla,Nle,Nt,T2_Cli,"t2_m")
      Deallocate (xlon,xlat,xled,xtime)
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,Pr_Cli_Obs,"pr_o")
      Deallocate (xlon,xlat,xtime)
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,T2_Cli_Obs,"t2_o")
      Deallocate (xlon,xlat,xtime)
   End Subroutine Read_Clim_Data
   !------------------------------------------------------------------------------
   Subroutine Read_Ope_Data
      Integer :: ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,icas, it1, it2
      Character (Len=255) :: Fname
      Real, Allocatable :: X4d(:,:,:,:), time(:), xlon(:),xlat(:),xmem(:),xtime(:)

      ! Pr_Ope(:,:,:,:,:,:), T2_Ope(:,:,:,:,:,:)     ! Nlon,Nlat,Nmem,Nlead,Ncase,Nmodel
      Allocate (T2_Ope(Nlon,Nlat,Nmem,Nlead,Ncase,Nmodel))
      Allocate (Pr_Ope(Nlon,Nlat,Nmem,Nlead,Ncase,Nmodel))
      Do imod=1,Nmodel
         write (Fname,"(A,I4.4,I2.2,A)") trim(Dir_Ope)//"Dr_Ope_",YOpe,MOpe,"_"//trim(TVar)//"_"//trim(Model_nm(imod))//".nc"
         Do icas=1,Ncase
            Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,trim(Cases(icas)))
            T2_Ope(:,:,:,:,icas,imod)=X4d(:,:,:,:)
            Deallocate (xlon,xlat,xmem,xtime,X4d)
         Enddo
      Enddo
      Do imod=1,Nmodel
         write (Fname,"(A,I4.4,I2.2,A)") trim(Dir_Ope)//"Dr_Ope_",YOpe,MOpe,"_"//trim(RVar)//"_"//trim(Model_nm(imod))//".nc"
         Do icas=1,Ncase
            Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,trim(Cases(icas)))
            Pr_Ope(:,:,:,:,icas,imod)=X4d(:,:,:,:)
            Deallocate (xlon,xlat,xmem,xtime,X4d)
         Enddo
      Enddo
   End Subroutine Read_Ope_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"
   !--------------------------------------------------------------------
   END ! Program