Program Dr_IDXs_SEACLID
  ! Calculating PDSI, SPI, SPEI 
  use netcdf

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

  Real, parameter    :: amssng = -99., Miss=amssng
  Integer, Parameter :: NVar = 2, Nmodels=7, Nseason=4, Nq=1000 
  Integer, Parameter :: NMonth = 12, Nmon=NMonth

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

  Real, Allocatable    :: beta_spi(:), gamm(:), pzero(:), SPI(:,:,:), SPE(:,:,:)
  Real, Allocatable, Dimension (:,:) :: mon_ori_pdsi, mon_ori_wplm, mon_ori_zind
  Real, Allocatable, Dimension (:,:) :: mon_sel_pdsi, mon_sel_wplm, mon_sel_zind

  Real, Allocatable :: Dat_Mon(:,:,:,:), DatNor(:,:,:)   ! For calculating PDSI
  Real, Allocatable :: AWC(:), xlon(:), xlat(:), Mask(:,:), xlon1(:), xlat1(:)
  Real, Allocatable :: Mask_Land(:), Lon(:), Lat(:)  ! (Nlon1*Nlat1) - All gridpoints
  !Real, Allocatable, Dimension (:,:) :: Prec, Temp
  Real, allocatable, dimension (:,:,:) :: Temp_Inp, Rain_Inp

  Real, Allocatable :: T2m(:,:,:), Pr(:,:,:)

!------------------------------------------------------------------------------
  Character (Len=255) :: Fname1,Fname, R_Fname, T_Fname, BL_Dir, PDSI_Dir, PDSI_WDir, Data_dir1,Data_dir2,Data_dir3, MaskFname, AWC_Fname, Out_Dir
  Character (Len=3)   :: TVar, RVar
  Character (Len=50)  :: Model_nm(Nmodels)
  Logical :: status, Empty

  Integer :: YBas1, YBas2, NYear, maxyrs
  Integer :: NLen, Nst, Nlon, Nlat, IO, Ntime, Nt, it
  Integer :: imod, ist, iyr, imw, ivar, isys, system, ilon, ilat
!------------------------------------------------------------------------------
  Real, Allocatable :: R_ORG(:,:,:,:),T_ORG(:,:,:,:),R_QM(:,:,:,:),T_QM(:,:,:,:)          ! (Nlon,Nlat,Nmodels,Ntime)
  Real, Allocatable :: R_ORG_S(:,:,:),T_ORG_S(:,:,:),R_QM_S(:,:,:),T_QM_S(:,:,:)  ! (Nmodels,Nseason,Landpoint*Ntime/4)

  Real, Allocatable :: PET_ORG_S(:,:,:),WB_ORG_S(:,:,:),PET_QM_S(:,:,:),WB_QM_S(:,:,:)  ! (Nmodels,Nseason,Landpoint*Ntime/4)
  Real, Allocatable :: PET_ORG_q(:,:,:),WB_ORG_q(:,:,:),PET_QM_q(:,:,:),WB_QM_q(:,:,:)  ! (Nmodels,Nseason,Nq)
  Real, Allocatable :: PET_ORG(:,:,:,:),WB_ORG(:,:,:,:),PET_QM(:,:,:,:),WB_QM(:,:,:,:)          ! (Nlon,Nlat,Nmodels,Ntime)

  Real, Allocatable, Dimension (:,:,:) :: Cor_R_ORG,Cor_T_ORG,Cor_R_QM,Cor_T_QM                     ! (Nlon,Nlat,Nmodels-1)
  Real, Allocatable, Dimension (:,:,:) :: prb_R_ORG,prb_T_ORG,prb_R_QM,prb_T_QM 
  Real, Allocatable, Dimension (:,:,:) :: Cor_WB_ORG,Cor_WB_QM
  Real, Allocatable, Dimension (:,:,:) :: prb_WB_ORG,prb_WB_QM

  Real, Allocatable, Dimension (:,:)   :: Ave_R_ORG,Ave_T_ORG,Ave_WB_ORG                            ! (Ntime,Nmodels)
  Real, Allocatable, Dimension (:,:)   :: Ave_R_QM, Ave_T_QM, Ave_WB_QM 

  Logical :: mms_to_monthly    ! True, False
  Character (Len=4) :: QMBC    ! "ORIG", "QM"
  Integer :: Landpoint,Ntime_S
!------------------------------------------------------------------------------
! For writing to NC file:
  Real, allocatable, Dimension (:, :, :) :: X1, X2, X3, X4, X5, X6
  Real, allocatable, Dimension (:) :: Xt
!------------------------------------------------------------------------------
  Namelist /File_Dir_Names/ Data_dir1,Data_dir2,Data_dir3, Model_nm, MaskFname, Out_Dir, &
                            TVar, RVar
  Namelist /Data_Info/      YBas1, YBas2
!------------------------------------------------------------------------------

  Open (1, file="Namelist_4season.txt", Status="Old")
  Read (1, File_Dir_Names)
  Read (1, Data_Info)
  
  Call Proc_land_mask_AWC    ! Return xlon, xlat, Nlon, Nlat, Nst, Mask, Mask_Land, Lon, Lat
  Nyear = (YBas2 - YBas1 + 1); Ntime = Nyear*NMonth; maxyrs = NYear
  Ntime_S=Landpoint*(Ntime/4)
  Call Memory_Allocate
  !print*,"Nlon,Nlat,Landpoint,Ntime,Ntime_S: ",Nlon,Nlat,Landpoint,Ntime,Ntime_S; pause

!  Allocate(Temp_Inp(1:Nlon,1:Nlat,1:Ntime));
!  Allocate(Rain_Inp(1:Nlon,1:Nlat,1:Ntime));

!  Allocate (R_ORG(Nlon,Nlat,Nmodels,Ntime),T_ORG(Nlon,Nlat,Nmodels,Ntime))
!  Allocate (R_QM(Nlon,Nlat,Nmodels,Ntime), T_QM(Nlon,Nlat,Nmodels,Ntime))

!  Allocate (PET_ORG_S(Nmodels,Nseason,Ntime_S),WB_ORG_S(Nmodels,Nseason,Ntime_S))
!  Allocate (PET_QM_S(Nmodels,Nseason,Ntime_S), WB_QM_S(Nmodels,Nseason,Ntime_S))
!  PET_ORG_S=Miss; WB_ORG_S=Miss; PET_QM_S=Miss; WB_QM_S=Miss

  QMBC="ORIG"; mms_to_monthly=.True.
  Do imod=1,Nmodels
     if (imod<Nmodels) then
        QMBC="ORIG"; mms_to_monthly=.True.
        R_Fname=trim(Data_dir1)//trim(RVar)//"_"//trim(Model_nm(imod))//"_"//trim(QMBC)//"_Hist.nc"
        T_Fname=trim(Data_dir1)//trim(TVar)//"_"//trim(Model_nm(imod))//"_"//trim(QMBC)//"_Hist.nc"
     else
        QMBC="APH"; mms_to_monthly=.False.
        R_Fname=trim(Data_dir3)//trim(RVar)//"_"//trim(Model_nm(imod))//"_Monthly_1976-2005.nc"
        T_Fname=trim(Data_dir3)//trim(TVar)//"_"//trim(Model_nm(imod))//"_Monthly_1976-2005.nc"
     endif
     Call Read_Pr_T2m
     print*,"Reading Data: Done: Nlat,Nlon,Nst=",Nlat,Nlon,Nst
     !
     Do ilat=1, Nlat
        Do ilon=1, Nlon
           Do it=1,Ntime
              R_ORG(ilon,ilat,imod,it)=Rain_Inp(ilon,ilat,it)
              T_ORG(ilon,ilat,imod,it)=Temp_Inp(ilon,ilat,it)
           Enddo
        Enddo
     Enddo
  Enddo
  QMBC="QM"; mms_to_monthly=.False.
  Do imod=1,Nmodels
     if (imod<Nmodels) then
        QMBC="QM"; mms_to_monthly=.False.
        R_Fname=trim(Data_dir2)//trim(RVar)//"_"//trim(Model_nm(imod))//"_"//trim(QMBC)//"_Hist.nc"
        T_Fname=trim(Data_dir2)//trim(TVar)//"_"//trim(Model_nm(imod))//"_"//trim(QMBC)//"_Hist.nc"
     else
        QMBC="APH"; mms_to_monthly=.False.
        R_Fname=trim(Data_dir3)//trim(RVar)//"_"//trim(Model_nm(imod))//"_Monthly_1976-2005.nc"
        T_Fname=trim(Data_dir3)//trim(TVar)//"_"//trim(Model_nm(imod))//"_Monthly_1976-2005.nc"
     endif
     Call Read_Pr_T2m
     print*,"Reading Data: Done: Nlat,Nlon,Nst=",Nlat,Nlon,Nst
     !
     Do ilat=1, Nlat
        Do ilon=1, Nlon
           Do it=1,Ntime
              R_QM(ilon,ilat,imod,it)=Rain_Inp(ilon,ilat,it)
              T_QM(ilon,ilat,imod,it)=Temp_Inp(ilon,ilat,it)
           Enddo
        Enddo
     Enddo
  Enddo
  Call Water_Balance
  Call Make_season_data
  Call Calc_Some_Stats
  Fname="NC1/Corr_Models_APH.nc"; Fname1="NC1/Areal_Mean_R_T_WB.nc"
  Call Write_Stat_to_NC (Fname,Fname1)
  print*,"Calling Calc_Q_Q "
  Call Calc_Q_Q
  Fname="NC1/PET_WB_season.nc"; Fname1="NC1/PET_WB_season_quantile.nc"
  Call Write_All_Dr_IDX_to_NetCDF(Fname,Fname1)
CONTAINS
!------------------------------------------------------------------------------
Subroutine Memory_Allocate
  Allocate(Temp_Inp(1:Nlon,1:Nlat,1:Ntime));
  Allocate(Rain_Inp(1:Nlon,1:Nlat,1:Ntime));

  Allocate (R_ORG(Nlon,Nlat,Nmodels,Ntime),T_ORG(Nlon,Nlat,Nmodels,Ntime))
  Allocate (R_QM(Nlon,Nlat,Nmodels,Ntime), T_QM(Nlon,Nlat,Nmodels,Ntime))

  Allocate (PET_ORG_S(Nmodels,Nseason,Ntime_S),WB_ORG_S(Nmodels,Nseason,Ntime_S))
  Allocate (PET_QM_S(Nmodels,Nseason,Ntime_S), WB_QM_S(Nmodels,Nseason,Ntime_S))
  PET_ORG_S=Miss; WB_ORG_S=Miss; PET_QM_S=Miss; WB_QM_S=Miss

   Allocate (PET_ORG_q(Nmodels,Nseason,Nq),WB_ORG_q(Nmodels,Nseason,Nq))
   Allocate (PET_QM_q(Nmodels,Nseason,Nq), WB_QM_q(Nmodels,Nseason,Nq))

   Allocate (PET_ORG(Nlon,Nlat,Nmodels,Ntime),WB_ORG(Nlon,Nlat,Nmodels,Ntime))
   Allocate (PET_QM(Nlon,Nlat,Nmodels,Ntime), WB_QM(Nlon,Nlat,Nmodels,Ntime))

   Allocate (Cor_R_ORG(Nlon,Nlat,Nmodels-1),Cor_T_ORG(Nlon,Nlat,Nmodels-1),Cor_R_QM(Nlon,Nlat,Nmodels-1),Cor_T_QM(Nlon,Nlat,Nmodels-1))
   Allocate (prb_R_ORG(Nlon,Nlat,Nmodels-1),prb_T_ORG(Nlon,Nlat,Nmodels-1),prb_R_QM(Nlon,Nlat,Nmodels-1),prb_T_QM(Nlon,Nlat,Nmodels-1))
   Allocate (Cor_WB_ORG(Nlon,Nlat,Nmodels-1),Cor_WB_QM(Nlon,Nlat,Nmodels-1))
   Allocate (prb_WB_ORG(Nlon,Nlat,Nmodels-1),prb_WB_QM(Nlon,Nlat,Nmodels-1))
   Allocate (Ave_R_ORG(Nmodels,Ntime),Ave_T_ORG(Nmodels,Ntime),Ave_WB_ORG(Nmodels,Ntime))
   Allocate (Ave_R_QM (Nmodels,Ntime),Ave_T_QM (Nmodels,Ntime),Ave_WB_QM (Nmodels,Ntime))

   Cor_R_ORG=Miss; Cor_T_ORG=Miss; Cor_R_QM=Miss; Cor_T_QM=Miss;
   prb_R_ORG=Miss; prb_T_ORG=Miss; prb_R_QM=Miss; prb_T_QM=Miss;
   Cor_WB_ORG=Miss; Cor_WB_ORG=Miss; Cor_WB_QM=Miss; Cor_T_QM=Miss;
   prb_WB_ORG=Miss; prb_WB_ORG=Miss; prb_WB_QM=Miss; prb_T_QM=Miss;
   Ave_R_ORG=Miss; Ave_T_ORG=Miss; Ave_WB_ORG=Miss;
   Ave_R_QM=Miss;  Ave_T_QM=Miss;  Ave_WB_QM=Miss;

End Subroutine Memory_Allocate
!------------------------------------------------------------------------------
Subroutine Calc_Some_Stats
   Integer :: ilon,ilat,imod,isea,ilan,it,it1,it2,it3,it4,imon
   Real, Allocatable :: X3d(:,:,:),Y3d(:,:,:),Co(:,:),Po(:,:)
   Real :: cor,prob

!   Allocate (Cor_R_ORG(Nlon,Nlat,Nmodels-1),Cor_T_ORG(Nlon,Nlat,Nmodels-1),Cor_R_QM(Nlon,Nlat,Nmodels-1),Cor_T_QM(Nlon,Nlat,Nmodels-1))
!   Allocate (prb_R_ORG(Nlon,Nlat,Nmodels-1),prb_T_ORG(Nlon,Nlat,Nmodels-1),prb_R_QM(Nlon,Nlat,Nmodels-1),prb_T_QM(Nlon,Nlat,Nmodels-1))
!   Allocate (Cor_WB_ORG(Nlon,Nlat,Nmodels-1),Cor_WB_ORG(Nlon,Nlat,Nmodels-1),Cor_WB_QM(Nlon,Nlat,Nmodels-1),Cor_T_QM(Nlon,Nlat,Nmodels-1))
!   Allocate (prb_WB_ORG(Nlon,Nlat,Nmodels-1),prb_WB_ORG(Nlon,Nlat,Nmodels-1),prb_WB_QM(Nlon,Nlat,Nmodels-1),prb_T_QM(Nlon,Nlat,Nmodels-1))
!   Allocate (Ave_R_ORG(Nmodels,Ntime),Ave_T_ORG(Nmodels,Ntime),Ave_WB_ORG(Nmodels,Ntime))
!   Cor_R_ORG=Miss; Cor_T_ORG=Miss; Cor_R_QM=Miss; Cor_T_QM=Miss; 
!   prb_R_ORG=Miss; prb_T_ORG=Miss; prb_R_QM=Miss; prb_T_QM=Miss; 
!   Cor_WB_ORG=Miss; Cor_WB_ORG=Miss; Cor_WB_QM=Miss; Cor_T_QM=Miss; 
!   prb_WB_ORG=Miss; prb_WB_ORG=Miss; prb_WB_QM=Miss; prb_T_QM=Miss; 
!   Ave_R_ORG=Miss; Ave_T_ORG=Miss; Ave_WB_ORG=Miss; 
!   Ave_R_QM=Miss;  Ave_T_QM=Miss;  Ave_WB_QM=Miss; 
!  Allocate (R_ORG(Nlon,Nlat,Nmodels,Ntime),T_ORG(Nlon,Nlat,Nmodels,Ntime))
!  Allocate (R_QM(Nlon,Nlat,Nmodels,Ntime), T_QM(Nlon,Nlat,Nmodels,Ntime))

   Allocate (X3d(Nlon,Nlat,Ntime),Y3d(Nlon,Nlat,Ntime))
   Allocate (Co(Nlon,Nlat),Po(Nlon,Nlat))
   Do imod=1,Nmodels-1
      Y3d(:,:,1:Ntime)=R_ORG(:,:,Nmodels,1:Ntime)
      X3d(:,:,1:Ntime)=R_ORG(:,:,imod,1:Ntime)
      Call Calc_Corr (X3d,Y3d,Co,Po); Cor_R_ORG(:,:,imod)=Co(:,:); prb_R_ORG(:,:,imod)=Po(:,:)

      Y3d(:,:,1:Ntime)=T_ORG(:,:,Nmodels,1:Ntime)
      X3d(:,:,1:Ntime)=T_ORG(:,:,imod,1:Ntime)
      Call Calc_Corr (X3d,Y3d,Co,Po); Cor_T_ORG(:,:,imod)=Co(:,:); prb_T_ORG(:,:,imod)=Po(:,:)

      Y3d(:,:,1:Ntime)=R_QM(:,:,Nmodels,1:Ntime)
      X3d(:,:,1:Ntime)=R_QM(:,:,imod,1:Ntime)
      Call Calc_Corr (X3d,Y3d,Co,Po); Cor_R_QM(:,:,imod)=Co(:,:); prb_R_QM(:,:,imod)=Po(:,:)

      Y3d(:,:,1:Ntime)=T_QM(:,:,Nmodels,1:Ntime)
      X3d(:,:,1:Ntime)=T_QM(:,:,imod,1:Ntime)
      Call Calc_Corr (X3d,Y3d,Co,Po); Cor_T_QM(:,:,imod)=Co(:,:); prb_T_QM(:,:,imod)=Po(:,:)

      Y3d(:,:,1:Ntime)=WB_ORG(:,:,Nmodels,1:Ntime)
      X3d(:,:,1:Ntime)=WB_ORG(:,:,imod,1:Ntime)
      Call Calc_Corr (X3d,Y3d,Co,Po); Cor_WB_ORG(:,:,imod)=Co(:,:); prb_WB_ORG(:,:,imod)=Po(:,:)

      Y3d(:,:,1:Ntime)=WB_QM(:,:,Nmodels,1:Ntime)
      X3d(:,:,1:Ntime)=WB_QM(:,:,imod,1:Ntime)
      Call Calc_Corr (X3d,Y3d,Co,Po); Cor_WB_QM(:,:,imod)=Co(:,:); prb_WB_QM(:,:,imod)=Po(:,:)
   Enddo
   Deallocate (X3d,Y3d,Po)
!   Allocate (Ave_R_ORG(Nmodels,Ntime),Ave_T_ORG(Nmodels,Ntime),Ave_WB_ORG(Nmodels,Ntime))
   Do imod=1,Nmodels
      Do it=1,Ntime
         Co(:,:)=R_ORG(:,:,imod,it); Ave_R_ORG(imod,it)=Areal_Ave(Co)
         Co(:,:)=T_ORG(:,:,imod,it); Ave_T_ORG(imod,it)=Areal_Ave(Co)
         Co(:,:)=R_QM (:,:,imod,it); Ave_R_QM (imod,it)=Areal_Ave(Co)
         Co(:,:)=T_QM (:,:,imod,it); Ave_T_QM (imod,it)=Areal_Ave(Co)
         Co(:,:)=WB_ORG(:,:,imod,it);Ave_WB_ORG(imod,it)=Areal_Ave(Co)
         Co(:,:)=WB_QM (:,:,imod,it);Ave_WB_QM (imod,it)=Areal_Ave(Co)
      Enddo
   Enddo
   Deallocate (Co)
End Subroutine Calc_Some_Stats
!------------------------------------------------------------------------------
Subroutine Write_Stat_to_NC (Fname1,Fname2)
   Character (Len=*) :: Fname1,Fname2
   Integer :: ist, jlat, ilon, it, ilen, ix, Nmod1
   Character (Len=250) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, t_nm
   Real :: xmodel(Nmodels)

   Nmod1=Nmodels-1
   if (Allocated(Xt)) Deallocate(Xt); Allocate (Xt(Nmod1))

   Do it=1,Nmod1;   Xt(it) = it; Enddo
   Write(t_unit,"(A,I4.4,A)") "months since ",YBas1,"-01-01"
   var_unit="-"; lon_nm="lon"; lat_nm="lat"; t_nm="Model"; Var_nm="r_org"
   Var_ln="Correlation between precipitation of Original model and APH. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI"
   Call wr_ncdf_3d (trim(Fname1),xlon,xlat,Xt,Nlon,Nlat,Nmod1,Cor_R_ORG,trim(Var_nm),t_unit,Var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
   Var_nm="p_r_org"
   Var_ln="P-Value of Correlation between precipitation of Original model and APH. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_3d_apd (trim(Fname1),prb_R_ORG,trim(Var_nm),Nlon,Nlat,Nmod1,Var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)

   Var_nm="r_qm"
   Var_ln="Correlation between precipitation of BC model and APH. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI"
   Call wr_ncdf_3d_apd (trim(Fname1),Cor_R_QM,trim(Var_nm),Nlon,Nlat,Nmod1,Var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
   Var_nm="p_r_qm"
   Var_ln="P-Value of Correlation between precipitation of BC model and APH. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_3d_apd (trim(Fname1),prb_R_QM,trim(Var_nm),Nlon,Nlat,Nmod1,Var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
   
   Var_nm="t_org"
   Var_ln="Correlation between temperature of Original model and APH. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI"
   Call wr_ncdf_3d_apd (trim(Fname1),Cor_T_ORG,trim(Var_nm),Nlon,Nlat,Nmod1,Var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
   Var_nm="p_t_org"
   Var_ln="P-Value of Correlation between temperature of Original model and APH. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_3d_apd (trim(Fname1),prb_T_ORG,trim(Var_nm),Nlon,Nlat,Nmod1,Var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)

   Var_nm="t_qm"
   Var_ln="Correlation between temperature of BC model and APH. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI"
   Call wr_ncdf_3d_apd (trim(Fname1),Cor_T_QM,trim(Var_nm),Nlon,Nlat,Nmod1,Var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
   Var_nm="p_t_qm"
   Var_ln="P-Value of Correlation between temperature of BC model and APH. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_3d_apd (trim(Fname1),prb_T_QM,trim(Var_nm),Nlon,Nlat,Nmod1,Var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)

   Var_nm="wb_org"
   Var_ln="Correlation between water balance of Original model and APH. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI"
   Call wr_ncdf_3d_apd (trim(Fname1),Cor_WB_ORG,trim(Var_nm),Nlon,Nlat,Nmod1,Var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
   Var_nm="p_wb_org"
   Var_ln="P-Value of Correlation between water balance of Original model and APH. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_3d_apd (trim(Fname1),prb_WB_ORG,trim(Var_nm),Nlon,Nlat,Nmod1,Var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)

   Var_nm="wb_qm"
   Var_ln="Correlation between water balance of BC model and APH. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI"
   Call wr_ncdf_3d_apd (trim(Fname1),Cor_WB_QM,trim(Var_nm),Nlon,Nlat,Nmod1,Var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
   Var_nm="p_wb_qm"
   Var_ln="P-Value of Correlation between water balance of BC model and APH. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_3d_apd (trim(Fname1),prb_WB_QM,trim(Var_nm),Nlon,Nlat,Nmod1,Var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)

   if (Allocated(Xt)) Deallocate(Xt); Allocate (Xt(Ntime))
   Do it=1,Ntime;   Xt(it) = it; Enddo
   Do it=1,Nmodels; xmodel(it)=it; Enddo
!   Allocate (Ave_R_ORG(Nmodels,Ntime),Ave_T_ORG(Nmodels,Ntime),Ave_WB_ORG(Nmodels,Ntime))

   lon_nm="model"; lat_nm="month"; 
   Var_nm="r_org"; var_unit="mm/month";
   var_ln="Areal mean of monthly original model precipitation. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_2d (Fname2,xmodel,Xt,Nmodels,Ntime,Ave_R_ORG,var_nm,var_ln,var_unit,lon_nm,lat_nm,Miss) 

   Var_nm="r_qm"; var_unit="mm/month";
   var_ln="Areal mean of monthly BC model precipitation. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_2d_apd (Fname2,Ave_R_QM,var_nm,Nmodels,Ntime,var_ln,var_unit,lon_nm,lat_nm,Miss)

   Var_nm="t_org"; var_unit="C";
   var_ln="Areal mean of monthly original model temperature. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_2d_apd (Fname2,Ave_T_ORG,var_nm,Nmodels,Ntime,var_ln,var_unit,lon_nm,lat_nm,Miss)

   Var_nm="t_qm"; var_unit="C";
   var_ln="Areal mean of monthly BC model temperature. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_2d_apd (Fname2,Ave_T_QM,var_nm,Nmodels,Ntime,var_ln,var_unit,lon_nm,lat_nm,Miss)

   Var_nm="wb_org"; var_unit="mm/month";
   var_ln="Areal mean of monthly original model water balance. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_2d_apd (Fname2,Ave_WB_ORG,var_nm,Nmodels,Ntime,var_ln,var_unit,lon_nm,lat_nm,Miss)

   Var_nm="wb_qm"; var_unit="mm/month";
   var_ln="Areal mean of monthly BC model water balance. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_2d_apd (Fname2,Ave_WB_QM,var_nm,Nmodels,Ntime,var_ln,var_unit,lon_nm,lat_nm,Miss)
End Subroutine Write_Stat_to_NC 
!------------------------------------------------------------------------------
Real Function Areal_Ave(X)
   Integer :: ilon,ilat,it
   Real :: X(Nlon,Nlat), tmp

   tmp=0.; it=0
   Do ilon=1,Nlon
      Do ilat=1,Nlat
         if (Mask(ilon,ilat)>0) then
            if (X(ilon,ilat)/=Miss) then
               tmp=tmp+X(ilon,ilat); it=it+1
            endif
         endif
      Enddo
   Enddo
   if (it>0) then; Areal_Ave=tmp/real(it); else; Areal_Ave=Miss;endif
End Function Areal_Ave
!------------------------------------------------------------------------------
Subroutine Calc_Corr (X,Y,Corre,Prob)
   Integer :: ilon,ilat,it
   Real :: corr,prb
   Real :: X1(Ntime),Y1(Ntime),X(Nlon,Nlat,Ntime),Y(Nlon,Nlat,Ntime)
   Real :: Corre(Nlon,Nlat),Prob(Nlon,Nlat)

   Corre=Miss; Prob=Miss
   Do ilon=1,Nlon
      Do ilat=1,Nlat
         if (Mask(ilon,ilat)>0) then
            X1(1:Ntime)=X(ilon,ilat,1:Ntime); Y1(1:Ntime)=Y(ilon,ilat,1:Ntime)
            Call Correlation (X1,Y1,Ntime,Corr,Prb,Miss)
            Corre(ilon,ilat)=Corr; Prob(ilon,ilat)=Prb
         endif
      Enddo
   Enddo
End Subroutine Calc_Corr 
!------------------------------------------------------------------------------
! thornthwaite (tempSeries, n, lat, etpSeries)
! Real, Allocatable :: PET_ORG(:,:,:,:),WB_ORG(:,:,:,:),PET_QM(:,:,:,:),WB_QM(:,:,:,:)          ! (Nlon,Nlat,Nmodels,Ntime)
Subroutine Water_Balance
   Integer :: ilon,ilat,imod,isea,ilan,it,it1,it2,it3,it4,imon
   Real, Allocatable :: Temp(:), ET(:)
   Real :: latx,petx

!   Allocate (PET_ORG(Nlon,Nlat,Nmodels,Ntime),WB_ORG(Nlon,Nlat,Nmodels,Ntime))
!   Allocate (PET_QM(Nlon,Nlat,Nmodels,Ntime), WB_QM(Nlon,Nlat,Nmodels,Ntime))
   Allocate (Temp(0:Ntime-1),ET(0:Ntime-1))
   PET_ORG=Miss; PET_QM=Miss; WB_ORG=Miss; WB_QM=Miss
   Do imod=1,Nmodels
      print*,"Water_Balance, imod: ",imod
      Do ilon=1,Nlon
         Do ilat=1,Nlat
            if (Mask(ilon,ilat)>0) then
               latx=xlat(ilat)
               Temp(0:Ntime-1)=T_ORG(ilon,ilat,imod,1:Ntime); Call thornthwaite(Temp,Ntime,latx,ET)
               PET_ORG(ilon,ilat,imod,1:Ntime)=ET(0:Ntime-1); 
               WB_ORG(ilon,ilat,imod,1:Ntime)=R_ORG(ilon,ilat,imod,1:Ntime)-PET_ORG(ilon,ilat,imod,1:Ntime)

               Temp(0:Ntime-1)=T_QM(ilon,ilat,imod,1:Ntime); Call thornthwaite(Temp,Ntime,latx,ET)
               PET_QM(ilon,ilat,imod,1:Ntime)=ET(0:Ntime-1); 
               WB_QM(ilon,ilat,imod,1:Ntime)=R_QM(ilon,ilat,imod,1:Ntime)-PET_QM(ilon,ilat,imod,1:Ntime)
            endif
         Enddo
      Enddo
   Enddo
End Subroutine Water_Balance
!------------------------------------------------------------------------------
Subroutine Make_season_data
   Integer :: ilon,ilat,imod,isea,ilan,it,it1,it2,it3,it4,imon

   print*,"Make_season_data..."
   Do imod=1,Nmodels
      print*,"Model: ",trim(Model_nm(imod))
      !Do isea=1,Nseason
         it1=0; it2=0; it3=0; it4=0;
         Do ilon=1,Nlon
            Do ilat=1,Nlat
               if (Mask(ilon,ilat)>0) then
                  Do it=1,Ntime
                     imon=mod(it,Nmon); if (imon==0) imon=Nmon
                     !print*,"imon: ", imon
                     Select Case (imon)
                        Case (12,1,2)
                           it1=it1+1;
                           PET_ORG_S(imod,1,it1)=PET_ORG(ilon,ilat,imod,it); WB_ORG_S(imod,1,it1)=WB_ORG(ilon,ilat,imod,it)
                           PET_QM_S (imod,1,it1)=PET_QM (ilon,ilat,imod,it); WB_QM_S (imod,1,it1)=WB_QM (ilon,ilat,imod,it)
                        Case (3,4,5)
                           it2=it2+1;
                           PET_ORG_S(imod,2,it2)=PET_ORG(ilon,ilat,imod,it); WB_ORG_S(imod,2,it2)=WB_ORG(ilon,ilat,imod,it)
                           PET_QM_S (imod,2,it2)=PET_QM (ilon,ilat,imod,it); WB_QM_S (imod,2,it2)=WB_QM (ilon,ilat,imod,it)
                        Case (6,7,8)
                           it3=it3+1;
                           PET_ORG_S(imod,3,it3)=PET_ORG(ilon,ilat,imod,it); WB_ORG_S(imod,3,it3)=WB_ORG(ilon,ilat,imod,it)
                           PET_QM_S (imod,3,it3)=PET_QM (ilon,ilat,imod,it); WB_QM_S (imod,3,it3)=WB_QM (ilon,ilat,imod,it)
                        Case (9,10,11)
                           it4=it4+1;
                           PET_ORG_S(imod,4,it4)=PET_ORG(ilon,ilat,imod,it); WB_ORG_S(imod,4,it4)=WB_ORG(ilon,ilat,imod,it)
                           PET_QM_S (imod,4,it4)=PET_QM (ilon,ilat,imod,it); WB_QM_S (imod,4,it4)=WB_QM (ilon,ilat,imod,it)
                     End Select 
                  Enddo
               endif
            Enddo
         Enddo
      !Enddo
      print*,it1,it2,it3,it4
   Enddo
!  Allocate (R_ORG_S(Nmodels,Nseason,Ntime_S),T_ORG_S(Nmodels,Nseason,Ntime_S))
!  Allocate (R_QM_S(Nmodels,Nseason,Ntime_S), T_QM_S(Nmodels,Nseason,Ntime_S))

End Subroutine Make_season_data
!------------------------------------------------------------------------------
Subroutine Proc_land_mask_AWC
   Integer :: ist, ilon, ilat, nlo, nla
   Real, Allocatable :: xtmp(:,:),lontmp(:),lattmp(:)

     Call rd_ncdf_2d (trim(MaskFname),xlon,xlat,Nlon,Nlat,Mask, "land")

     Nst=Nlon*Nlat
     Allocate (Mask_Land(Nst))
     Allocate (lon(Nst), lat(Nst)) !, AWC(Nst))
     ist=0; Landpoint=0
     Do ilat=1, Nlat
        Do ilon=1, Nlon
           ist=ist+1; Mask_Land(ist)=Mask(ilon,ilat);
           if (Mask(ilon,ilat)>0) Landpoint=Landpoint+1
!           AWC(ist)=xtmp(ilon,ilat)
           Lon(ist)=xlon(ilon); Lat(ist)=xlat(ilat)
        Enddo
     Enddo
!     Deallocate (lontmp,lattmp,xtmp) 
End Subroutine Proc_land_mask_AWC
!------------------------------------------------------------------------------
Subroutine Read_Pr_T2m
  Real, Allocatable :: XTime(:)
  Real, Allocatable :: xlon(:), xlat(:)
  Integer :: Nlon1, Nlat1

  print*,"==================>In Read_Pr_T2m: "
  print*,"Reading ",trim(R_Fname)," and ", trim(T_Fname) 
  Call rd_ncdf_3d (trim(R_Fname), xlon, xlat, XTime, Nlon1, Nlat1, Nt, Pr, trim(RVar))
  Deallocate (xlon, xlat, XTime)
  Call rd_ncdf_3d (trim(T_Fname), xlon, xlat, XTime, Nlon1, Nlat1, Nt, T2m, trim(TVar))
  print*,"Nlon1, Nlat1, Nt=",Nlon1, Nlat1, Nt
  Call Proc_Pr_T2m 
  Deallocate (Pr, T2m)
End Subroutine Read_Pr_T2m
!------------------------------------------------------------------------------
Subroutine Proc_Pr_T2m 
  Integer, parameter :: NMon_Skip = 0 !72
  Real :: TMean(Nmonth),RMean(Nmonth),xtas(Nyear),xpr(Nyear),xt75(Nmonth),xp75(Nmonth)
  Integer :: iyr,i,j,it,ict,icr,imon,imon1
  Real, Allocatable :: Xtmp(:,:)
  
  print*,"Beginning Proc_Pr_T2m...:", Ntime,Nt,NMon_Skip 
  Where(T2m==-99.9) T2m=Miss; Temp_Inp(:,:,1:Ntime)=T2m(:,:,NMon_Skip+1:Nt); ! SL tu 1976 - 2005
  Where(Pr ==-99.9) Pr =Miss; Rain_Inp(:,:,1:Ntime)=Pr(:,:,NMon_Skip+1:Nt)  ! Monthly data, mm/s

  ! Fill up missing data:
  Do j=1, Nlat
     Do i=1, Nlon
        if (Mask(i,j) == 1) then
           Do imon=1, Nmonth
              TMean(imon)=0; RMean(imon)=0; ict=0; icr=0
              Do it=1, Ntime
                 imon1=mod(it,Nmonth); if (imon1==0) imon1=Nmonth
                 if (imon1==imon) then
                    if (Temp_Inp(i,j,it)/=Miss) then
                       !TMean(imon)=TMean(imon)+Temp_Inp(i,j,it); ict=ict+1
                       ict=ict+1; xtas(ict)=Temp_Inp(i,j,it)
                    endif
                    if (Rain_Inp(i,j,it)>=0) then !/=Miss) then
                       !RMean(imon)=RMean(imon)+Rain_Inp(i,j,it); icr=icr+1
                       icr=icr+1; xpr(icr)=Rain_Inp(i,j,it)
                    endif
                 endif
              Enddo
              !TMean(imon)=TMean(imon)/Real(ict)
              !RMean(imon)=RMean(imon)/Real(icr)
              if (ict>Nyear/3*2.AND.icr>Nyear/3*2) then
                 TMean(imon)=Percentile(xtas,ict,0.5,Miss)
                 RMean(imon)=Percentile(xpr, icr,0.5,Miss)
                 xt75 (imon)=Percentile(xtas,ict,0.75,Miss)
                 xp75 (imon)=Percentile(xpr, icr,0.75,Miss)
              else
                 TMean(imon)=Miss; RMean(imon)=Miss; xt75 (imon)=Miss; xp75 (imon)=Miss
              endif
           Enddo
           Do it=1, Ntime
              imon=mod(it,Nmonth); if (imon==0) imon=Nmonth
              if (Temp_Inp(i,j,it)==Miss) Temp_Inp(i,j,it)=TMean(imon)
              if (Rain_Inp(i,j,it)==Miss) Rain_Inp(i,j,it)=RMean(imon)
              if (Temp_Inp(i,j,it)>2.0*xt75(imon)) Temp_Inp(i,j,it)=xt75(imon)
              if (Rain_Inp(i,j,it)>3.0*xp75(imon)) Rain_Inp(i,j,it)=xp75(imon)
           Enddo
        endif
     Enddo
  Enddo
  ! Convert Pr from mm/s to mm/month and K to C:
  Allocate (Xtmp(1:Nlon,1:Nlat))
  it=0
  Do iyr=YBas1, YBas2
     Do imon=1, Nmonth
        it=it+1
        Xtmp(:,:)=Rain_Inp(:,:,it)
        !Where (Xtmp/=Miss) Xtmp=Xtmp*DayOfMon(imon,iyr)  ! Daily to Monthly
        if (mms_to_monthly) then
           Where (Xtmp/=Miss) Xtmp=Xtmp*86400.0*DayOfMon(imon,iyr)  ! mm/s to Monthly
        endif
        Where (Xtmp/=Miss) Xtmp=Xtmp;
        Rain_Inp(:,:,it)=Xtmp(:,:)
        Xtmp(:,:)=Temp_Inp(:,:,it)
        if (mms_to_monthly) then
           Where (Xtmp/=Miss) Xtmp=Xtmp-273.15    ! K to C
        endif
        Temp_Inp(:,:,it)=Xtmp(:,:)
     Enddo
  Enddo
  ! Assign Miss values to all Sea-points:
  Do it=1, Ntime
     Do j=1, Nlat
        Do i=1, Nlon
           if (Mask(i,j)==Miss) then
              Rain_Inp(i,j,it)=Miss; Temp_Inp(i,j,it)=Miss
           endif
        Enddo
     Enddo
  Enddo
  Deallocate(Xtmp)
End Subroutine Proc_Pr_T2m
!------------------------------------------------------------------------------
Subroutine Calc_Q_Q
   Integer :: imod,isea,it,i,j,k
   Real :: q(Nq), p
   Real, Allocatable :: X1d(:)

!  Real, Allocatable :: R_ORG_q(:,:,:),T_ORG_q(:,:,:),R_QM_q(:,:,:),T_QM_q(:,:,:)  ! (Nmodels,Nseason,Nq)
! x = percentile(Y,N1,p,Miss)

   !print*,"In Calc_Q_Q....... ",Ntime_S; pause
   Allocate (X1d(Ntime_S))
!   Allocate (PET_ORG_q(Nmodels,Nseason,Nq),WB_ORG_q(Nmodels,Nseason,Nq))
!   Allocate (PET_QM_q(Nmodels,Nseason,Nq), WB_QM_q(Nmodels,Nseason,Nq))

   Do i=1,Nq; q(i)=Real(i)/Real(Nq+1); Enddo
   Do imod=1,Nmodels
      Do isea=1,Nseason
         print*,imod,isea
         print*,"1"
         X1d(1:Ntime_S)=PET_ORG_S(imod,isea,1:Ntime_S)
         Do i=1,Nq
            p = q(i)
!print*,p,Ntime_S,X1d(100:110)
            PET_ORG_q(imod,isea,i) = percentile(X1d,Ntime_S,p,Miss)
         Enddo
         print*,"2"
         X1d(1:Ntime_S)=WB_ORG_S(imod,isea,1:Ntime_S)
         Do i=1,Nq
            p = q(i)
            WB_ORG_q(imod,isea,i) = percentile(X1d,Ntime_S,p,Miss)
         Enddo
         print*,"3"
         X1d(1:Ntime_S)=PET_QM_S(imod,isea,1:Ntime_S)
         Do i=1,Nq
            p = q(i)
            PET_QM_q(imod,isea,i) = percentile(X1d,Ntime_S,p,Miss)
         Enddo
         print*,"4"
         X1d(1:Ntime_S)=WB_QM_S(imod,isea,1:Ntime_S)
         Do i=1,Nq
            p = q(i)
            WB_QM_q(imod,isea,i) = percentile(X1d,Ntime_S,p,Miss)
         Enddo
      Enddo
   Enddo
End Subroutine Calc_Q_Q
!------------------------------------------------------------------------------
Subroutine Write_All_Dr_IDX_to_NetCDF(Fname, Fname1) !, Nlon, Nlat, Ntime)
   Integer :: ist, jlat, ilon, it, ilen, ix
   Character (Len=250) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, t_nm
   Character (Len=*) :: Fname, Fname1
   Real :: X_season(Nseason), X_model(Nmodels)

!  Allocate (R_ORG_S(Nmodels,Nseason,Ntime_S),T_ORG_S(Nmodels,Nseason,Ntime_S))
!  Allocate (R_QM_S(Nmodels,Nseason,Ntime_S), T_QM_S(Nmodels,Nseason,Ntime_S))

   !if (.Not.Allocated(X1)) Allocate(X1(Nlon, Nlat, Ntime))
   !if (.Not.Allocated(X2)) Allocate(X2(Nlon, Nlat, Ntime))
!print*,R_ORG_S(:,:,100:110)
!print*,R_QM_S(:,:,100:110)
!print*,T_ORG_S(:,:,100:110)
!print*,T_QM_S(:,:,100:110)

   if (.Not.Allocated(Xt)) Allocate(Xt(Ntime_S))

   Do it=1, Ntime_S;   Xt(it) = it; Enddo
   Do it=1, Nseason; X_season(it) = it; Enddo
   Do it=1, Nmodels; X_model(it) = it; Enddo

   Write(t_unit,"(A,I4.4,A)") "months since ",YBas1,"-01-01"
   var_unit="mm/month"; lon_nm="Model"; lat_nm="Season"; t_nm="month"; Var_nm="e_org"
   Var_ln="Original model PET in the seasons. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_3d (trim(Fname),X_model,X_season,Xt,Nmodels,Nseason,Ntime_S,PET_ORG_S,trim(Var_nm),trim(t_unit),trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)
   Var_nm="e_qm"
   Var_ln="QM BC model PET in the seasons. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_3d_apd (trim(Fname),PET_QM_S,trim(Var_nm),Nmodels,Nseason,Ntime_S,trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)
   Var_nm="wb_org"
   Var_ln="Original model water balance in the seasons. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_3d_apd (trim(Fname),WB_ORG_S,trim(Var_nm),Nmodels,Nseason,Ntime_S,trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)
   Var_nm="wb_qm"
   Var_ln="QM BC model water balance in the seasons. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_3d_apd (trim(Fname),WB_QM_S,trim(Var_nm),Nmodels,Nseason,Ntime_S,trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)

   Deallocate (Xt); Allocate(Xt(Nq))
   Do it=1, Nq;   Xt(it) = it; Enddo

   Write(t_unit,"(A,I4.4,A)") "months since ",YBas1,"-01-01"
   var_unit="mm/month"; lon_nm="Model"; lat_nm="Season"; t_nm="quantile"; Var_nm="e_org"
   Var_ln="Original model PET in the seasons. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_3d (trim(Fname1),X_model,X_season,Xt,Nmodels,Nseason,Nq,PET_ORG_q,trim(Var_nm),trim(t_unit),trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)
   Var_nm="e_qm"
   Var_ln="QM BC model PET in the seasons. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_3d_apd (trim(Fname1),PET_QM_q,trim(Var_nm),Nmodels,Nseason,Nq,trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)
   Var_nm="wb_org"
   Var_ln="Original model water balance in the seasons. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_3d_apd (trim(Fname1),WB_ORG_q,trim(Var_nm),Nmodels,Nseason,Nq,trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)
   Var_nm="wb_qm"
   Var_ln="QM BC model water balance in the seasons. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH"
   Call wr_ncdf_3d_apd (trim(Fname1),WB_QM_q,trim(Var_nm),Nmodels,Nseason,Nq,trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng)

   print*,"************ Write to file Done !"
End Subroutine Write_All_Dr_IDX_to_NetCDF
!------------------------------------------------------------------------------
  Integer Function DayOfMon(Mon, Year)
  Integer :: Mon, Year
  Select Case (Mon)
     Case (1,3,5,7,8,10,12)
        DayOfMon = 31
     Case (4,6,9,11)
        DayOfMon = 30
     Case (2)
        if (mod(Year,4)==0.AND.mod(Year,100)/=0) then
           DayOfMon = 29
        else
           DayOfMon = 28
        endif
  End Select
  End Function DayOfMon
!--------------------------------------------------------------
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 "Percentile.f90"
Include "/home/tanpv/Mylibs/Calc_Correlation_Coefficient.f90"
!--------------------------------------------------------------------
END ! Program
Include "/home/tanpv/Mylibs/Student_Prob.f90"