Subroutine Calc_Drought_Characteristics

Integer :: ilon,ilat,ifile,imod,isce,it,ireg, imon,imon1, if1
Real :: ND1, ND2, ND3
Real :: X(Nmon)

! Sum over Ntime of the Number of Drought Month D1,2,3
Do ifile=1, Nfile
   Do ilon=1, Nlon
      Do ilat=1, Nlat
         if (Reg(ilon,ilat)>=1.AND.Reg(ilon,ilat)<=Nreg) then
            ND1=0; ND2=0; ND3=0
            Do it=1, Ntime
               if (D1_All(ilon,ilat,it,ifile)==1.) ND1=ND1+1
               if (D2_All(ilon,ilat,it,ifile)==1.) ND2=ND2+1
               if (D3_All(ilon,ilat,it,ifile)==1.) ND3=ND3+1
            Enddo
            Nmon_D1(ilon,ilat,ifile)=ND1; Nmon_D2(ilon,ilat,ifile)=ND2; Nmon_D3(ilon,ilat,ifile)=ND3 
         endif
      Enddo
   Enddo
Enddo
! Ensamble of the Number of Drought Month:
  Call Ensemble_1 (Nmon_D1); Call Ensemble_1 (Nmon_D2); Call Ensemble_1 (Nmon_D3);

! Sum over Ntime of Drought Severity
Do ifile=1, Nfile
   Do ilon=1, Nlon
      Do ilat=1, Nlat
         if (Reg(ilon,ilat)>=1.AND.Reg(ilon,ilat)<=Nreg) then
            ND1=0; ND2=0; ND3=0
            Do it=1, Ntime
            !   if (D1_All(ilon,ilat,it,ifile)==1.) ND1=ND1+PDSI_All(ilon,ilat,it,ifile)
            !   if (D2_All(ilon,ilat,it,ifile)==1.) ND2=ND2+PDSI_All(ilon,ilat,it,ifile)
            !   if (D3_All(ilon,ilat,it,ifile)==1.) ND3=ND3+PDSI_All(ilon,ilat,it,ifile)
               if (D1_All(ilon,ilat,it,ifile)==1. .AND. PDSI_All(ilon,ilat,it,ifile)/=Miss) ND1=ND1+PDSI_All(ilon,ilat,it,ifile)
               if (D2_All(ilon,ilat,it,ifile)==1. .AND. PDSI_All(ilon,ilat,it,ifile)/=Miss) ND2=ND2+PDSI_All(ilon,ilat,it,ifile)
               if (D3_All(ilon,ilat,it,ifile)==1. .AND. PDSI_All(ilon,ilat,it,ifile)/=Miss) ND3=ND3+PDSI_All(ilon,ilat,it,ifile) 
            Enddo
            Sum_DS1(ilon,ilat,ifile)=ND1; Sum_DS2(ilon,ilat,ifile)=ND2; Sum_DS3(ilon,ilat,ifile)=ND3
         endif
      Enddo
   Enddo
Enddo
! Ensamble of Drought Severity:
  Call Ensemble_1 (Sum_DS1); Call Ensemble_1 (Sum_DS2); Call Ensemble_1 (Sum_DS3);

! Monthly mean of Geo_Ext = Sum over Ntime of Geo_D1,2,3 divide to Ntime
Do ifile=1, Nfile
   Do ireg=1, Nreg+1
      Do imon=1, Nmon
         ND1=0; ND2=0; ND3=0
         Do it=1, Ntime
            imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon
            if (imon1==imon) then
               ND1=ND1+Geo_D1(ireg,it,ifile); ND2=ND2+Geo_D2(ireg,it,ifile); ND3=ND3+Geo_D3(ireg,it,ifile);
            endif
         Enddo
         Sum_Geo_D1(ireg,imon,ifile)=ND1/Real(NYear); Sum_Geo_D2(ireg,imon,ifile)=ND2/Real(NYear); Sum_Geo_D3(ireg,imon,ifile)=ND3/Real(NYear)
      Enddo
   Enddo
Enddo
! Ensemble of Geo_Ext
  Call Ensemble_2 (Sum_Geo_D1); Call Ensemble_2 (Sum_Geo_D2); Call Ensemble_2 (Sum_Geo_D3);

! Ensamble of Drought duration & inter-arrival time:
  Call Ensemble_1 (Mean_Duration_D1);  Call Ensemble_1 (Mean_Duration_D2);  Call Ensemble_1 (Mean_Duration_D3);
  Call Ensemble_1 (Inter_Arr_Time_D1); Call Ensemble_1 (Inter_Arr_Time_D2); Call Ensemble_1 (Inter_Arr_Time_D3);

! Fig.8,9: Changes of Number of Drought months between Scenerios and BL
! Fig. 10,11: Changes of Drought Severity between Scenerios and BL
! Fig. 12,13: Changes in mean inter-arrival time  
Do ilon=1, Nlon
   Do ilat=1, Nlat
   if (Reg(ilon,ilat)>0) then
      Do isce=1, NSCs
         Do imod=1, Nmodels+1
            ifile=isce*(Nmodels+1) + imod + 1
          ! Fig.8,9: Changes of Number of Drought months
          !  DelNmon_D1(ilon,ilat,imod,isce) = (Nmon_D1(ilon,ilat,ifile)-Nmon_D1(ilon,ilat,imod+1))/Nmon_D1(ilon,ilat,imod+1)*100.
          !  DelNmon_D2(ilon,ilat,imod,isce) = (Nmon_D2(ilon,ilat,ifile)-Nmon_D2(ilon,ilat,imod+1))/Nmon_D2(ilon,ilat,imod+1)*100.
          !  DelNmon_D3(ilon,ilat,imod,isce) = (Nmon_D3(ilon,ilat,ifile)-Nmon_D3(ilon,ilat,imod+1))/Nmon_D3(ilon,ilat,imod+1)*100.

            DelNmon_D1(ilon,ilat,imod,isce) = Changes (Nmon_D1(ilon,ilat,ifile), Nmon_D1(ilon,ilat,imod+1))
            DelNmon_D2(ilon,ilat,imod,isce) = Changes (Nmon_D2(ilon,ilat,ifile), Nmon_D2(ilon,ilat,imod+1))
            DelNmon_D3(ilon,ilat,imod,isce) = Changes (Nmon_D3(ilon,ilat,ifile), Nmon_D3(ilon,ilat,imod+1))

          ! Fig. 10,11: Changes of Drought Severity  
          !  DelSum_DS1(ilon,ilat,imod,isce) = (Sum_DS1(ilon,ilat,ifile)-Sum_DS1(ilon,ilat,imod+1))/Sum_DS1(ilon,ilat,imod+1)*100.
          !  DelSum_DS2(ilon,ilat,imod,isce) = (Sum_DS2(ilon,ilat,ifile)-Sum_DS2(ilon,ilat,imod+1))/Sum_DS2(ilon,ilat,imod+1)*100.
          !  DelSum_DS3(ilon,ilat,imod,isce) = (Sum_DS3(ilon,ilat,ifile)-Sum_DS3(ilon,ilat,imod+1))/Sum_DS3(ilon,ilat,imod+1)*100.

            DelSum_DS1(ilon,ilat,imod,isce) = Changes (Sum_DS1(ilon,ilat,ifile), Sum_DS1(ilon,ilat,imod+1))
            DelSum_DS2(ilon,ilat,imod,isce) = Changes (Sum_DS2(ilon,ilat,ifile), Sum_DS2(ilon,ilat,imod+1))
            DelSum_DS3(ilon,ilat,imod,isce) = Changes (Sum_DS3(ilon,ilat,ifile), Sum_DS3(ilon,ilat,imod+1))

          ! Fig. 12,13: Changes in mean inter-arrival time
          !  Del_Inter_Arr_Time_D1(ilon,ilat,imod,isce) = (Inter_Arr_Time_D1(ilon,ilat,ifile)-Inter_Arr_Time_D1(ilon,ilat,imod+1))/Inter_Arr_Time_D1(ilon,ilat,imod+1)*100.
          !  Del_Inter_Arr_Time_D2(ilon,ilat,imod,isce) = (Inter_Arr_Time_D2(ilon,ilat,ifile)-Inter_Arr_Time_D2(ilon,ilat,imod+1))/Inter_Arr_Time_D2(ilon,ilat,imod+1)*100.
          !  Del_Inter_Arr_Time_D3(ilon,ilat,imod,isce) = (Inter_Arr_Time_D3(ilon,ilat,ifile)-Inter_Arr_Time_D3(ilon,ilat,imod+1))/Inter_Arr_Time_D3(ilon,ilat,imod+1)*100.

            Del_Inter_Arr_Time_D1(ilon,ilat,imod,isce) = Changes (Inter_Arr_Time_D1(ilon,ilat,ifile), Inter_Arr_Time_D1(ilon,ilat,imod+1))
            Del_Inter_Arr_Time_D2(ilon,ilat,imod,isce) = Changes (Inter_Arr_Time_D2(ilon,ilat,ifile), Inter_Arr_Time_D2(ilon,ilat,imod+1))
            Del_Inter_Arr_Time_D3(ilon,ilat,imod,isce) = Changes (Inter_Arr_Time_D3(ilon,ilat,ifile), Inter_Arr_Time_D3(ilon,ilat,imod+1))

          ! Fig. 12b,13b: Changes in mean duration
          !  Del_Mean_Duration_D1(ilon,ilat,imod,isce) = (Mean_Duration_D1(ilon,ilat,ifile)-Mean_Duration_D1(ilon,ilat,imod+1))/Mean_Duration_D1(ilon,ilat,imod+1)*100.
          !  Del_Mean_Duration_D2(ilon,ilat,imod,isce) = (Mean_Duration_D2(ilon,ilat,ifile)-Mean_Duration_D2(ilon,ilat,imod+1))/Mean_Duration_D2(ilon,ilat,imod+1)*100.
          !  Del_Mean_Duration_D3(ilon,ilat,imod,isce) = (Mean_Duration_D3(ilon,ilat,ifile)-Mean_Duration_D3(ilon,ilat,imod+1))/Mean_Duration_D3(ilon,ilat,imod+1)*100.

            Del_Mean_Duration_D1(ilon,ilat,imod,isce) = Changes (Mean_Duration_D1(ilon,ilat,ifile), Mean_Duration_D1(ilon,ilat,imod+1))
            Del_Mean_Duration_D2(ilon,ilat,imod,isce) = Changes (Mean_Duration_D2(ilon,ilat,ifile), Mean_Duration_D2(ilon,ilat,imod+1))
            Del_Mean_Duration_D3(ilon,ilat,imod,isce) = Changes (Mean_Duration_D3(ilon,ilat,ifile), Mean_Duration_D3(ilon,ilat,imod+1))

          ! Tan: Changes of Number of Drought Spells
!            DelNSpel_D1(ilon,ilat,imod,isce) = (N_Spel_D1(ilon,ilat,ifile)-N_Spel_D1(ilon,ilat,imod+1))/N_Spel_D1(ilon,ilat,imod+1)*100.
!            DelNSpel_D2(ilon,ilat,imod,isce) = (N_Spel_D2(ilon,ilat,ifile)-N_Spel_D2(ilon,ilat,imod+1))/N_Spel_D2(ilon,ilat,imod+1)*100.
!            DelNSpel_D3(ilon,ilat,imod,isce) = (N_Spel_D3(ilon,ilat,ifile)-N_Spel_D3(ilon,ilat,imod+1))/N_Spel_D3(ilon,ilat,imod+1)*100.

          !  DelNSpel_D1(ilon,ilat,imod,isce) = (NSpell_D1(ilon,ilat,ifile)-NSpell_D1(ilon,ilat,imod+1))/NSpell_D1(ilon,ilat,imod+1)*100.
          !  DelNSpel_D2(ilon,ilat,imod,isce) = (NSpell_D2(ilon,ilat,ifile)-NSpell_D2(ilon,ilat,imod+1))/NSpell_D2(ilon,ilat,imod+1)*100.
          !  DelNSpel_D3(ilon,ilat,imod,isce) = (NSpell_D3(ilon,ilat,ifile)-NSpell_D3(ilon,ilat,imod+1))/NSpell_D3(ilon,ilat,imod+1)*100.

            DelNSpel_D1(ilon,ilat,imod,isce) = Changes (NSpell_D1(ilon,ilat,ifile), NSpell_D1(ilon,ilat,imod+1))
            DelNSpel_D2(ilon,ilat,imod,isce) = Changes (NSpell_D2(ilon,ilat,ifile), NSpell_D2(ilon,ilat,imod+1))
            DelNSpel_D3(ilon,ilat,imod,isce) = Changes (NSpell_D3(ilon,ilat,ifile), NSpell_D3(ilon,ilat,imod+1))
         Enddo
      Enddo
   endif
   Enddo
Enddo

! Fig. 14: Change in Spatial extent between Scenerios and BL
Do ireg=1, Nreg+1
   Do imon=1, Nmon
      Do isce=1, NSCs
         Do imod=1, Nmodels+1     ! ifile=1: OBS; ifile=2..Nmodels+2: BL; ifile=Nmodels+3..Nfile: Scenarios
            ifile=isce*(Nmodels+1) + imod + 1
            ! DelSum_Geo_D1(ireg,imon,imod,isce) = (Sum_Geo_D1(ireg,imon,ifile)-Sum_Geo_D1(ireg,imon,imod+1))/Sum_Geo_D1(ireg,imon,imod+1)*100. ! imod+1=2..8: BL 
            ! DelSum_Geo_D2(ireg,imon,imod,isce) = (Sum_Geo_D2(ireg,imon,ifile)-Sum_Geo_D2(ireg,imon,imod+1))/Sum_Geo_D2(ireg,imon,imod+1)*100. 
            ! DelSum_Geo_D3(ireg,imon,imod,isce) = (Sum_Geo_D3(ireg,imon,ifile)-Sum_Geo_D3(ireg,imon,imod+1))/Sum_Geo_D3(ireg,imon,imod+1)*100. 
            DelSum_Geo_D1(ireg,imon,imod,isce) = Changes (Sum_Geo_D1(ireg,imon,ifile), Sum_Geo_D1(ireg,imon,imod+1))
            DelSum_Geo_D2(ireg,imon,imod,isce) = Changes (Sum_Geo_D2(ireg,imon,ifile), Sum_Geo_D2(ireg,imon,imod+1))
            DelSum_Geo_D3(ireg,imon,imod,isce) = Changes (Sum_Geo_D3(ireg,imon,ifile), Sum_Geo_D3(ireg,imon,imod+1))
         Enddo
      Enddo
   Enddo
Enddo

! Regional mean of Drought Characteristics over entire periods and their CHANGES:

Call Regional_Mean (Nmon_D1, DMon1);           Call Regional_Mean (Nmon_D2, DMon2);           Call Regional_Mean (Nmon_D3, DMon3);
Call Regional_Mean (Mean_Duration_D1, DDur1);  Call Regional_Mean (Mean_Duration_D2, DDur2);  Call Regional_Mean (Mean_Duration_D3, DDur3)
Call Regional_Mean (Inter_Arr_Time_D1, DInt1); Call Regional_Mean (Inter_Arr_Time_D2, DInt2); Call Regional_Mean (Inter_Arr_Time_D3, DInt3)
Call Regional_Mean (Sum_DS1, DSev1);           Call Regional_Mean (Sum_DS2, DSev2);           Call Regional_Mean (Sum_DS3, DSev3)
!Call Regional_Mean (N_Spel_D1, DSpl1);         Call Regional_Mean (N_Spel_D2, DSpl2);         Call Regional_Mean (N_Spel_D3, DSpl3)
Call Regional_Mean (NSpell_D1, DSpl1);         Call Regional_Mean (NSpell_D2, DSpl2);         Call Regional_Mean (NSpell_D3, DSpl3)

Do ifile=1, Nfile
   Do ireg=1, Nreg+1
      X(1:Nmon)=Sum_Geo_D1(ireg,1:Nmon,ifile); DGeo1(ireg,ifile)=Time_Ave(X, Nmon)
      X(1:Nmon)=Sum_Geo_D2(ireg,1:Nmon,ifile); DGeo2(ireg,ifile)=Time_Ave(X, Nmon)
      X(1:Nmon)=Sum_Geo_D3(ireg,1:Nmon,ifile); DGeo3(ireg,ifile)=Time_Ave(X, Nmon)
   Enddo
Enddo

! Fig. 15->??: Change in Drought Characteristics between Scenerios and BL
Do ireg=1, Nreg+1
   Do isce=1, NSCs
      Do imod=1, Nmodels+1     ! ifile=1: OBS; ifile=2..Nmodels+2: BL; ifile=Nmodels+3..Nfile: Scenarios
         ifile=isce*(Nmodels+1) + imod + 1
         ! Del_DMon1(ireg,imod,isce) = (DMon1(ireg,ifile)-DMon1(ireg,imod+1))/DMon1(ireg,imod+1)*100.
         ! Del_DMon2(ireg,imod,isce) = (DMon2(ireg,ifile)-DMon2(ireg,imod+1))/DMon2(ireg,imod+1)*100.
         ! Del_DMon3(ireg,imod,isce) = (DMon3(ireg,ifile)-DMon3(ireg,imod+1))/DMon3(ireg,imod+1)*100.

         Del_DMon1(ireg,imod,isce) = Changes (DMon1(ireg,ifile), DMon1(ireg,imod+1))
         Del_DMon2(ireg,imod,isce) = Changes (DMon2(ireg,ifile), DMon2(ireg,imod+1))
         Del_DMon3(ireg,imod,isce) = Changes (DMon3(ireg,ifile), DMon3(ireg,imod+1))

         ! Del_DGeo1(ireg,imod,isce) = (DGeo1(ireg,ifile)-DGeo1(ireg,imod+1))/DGeo1(ireg,imod+1)*100.
         ! Del_DGeo2(ireg,imod,isce) = (DGeo2(ireg,ifile)-DGeo2(ireg,imod+1))/DGeo2(ireg,imod+1)*100.
         ! Del_DGeo3(ireg,imod,isce) = (DGeo3(ireg,ifile)-DGeo3(ireg,imod+1))/DGeo3(ireg,imod+1)*100.

         Del_DGeo1(ireg,imod,isce) = Changes (DGeo1(ireg,ifile), DGeo1(ireg,imod+1))
         Del_DGeo2(ireg,imod,isce) = Changes (DGeo2(ireg,ifile), DGeo2(ireg,imod+1))
         Del_DGeo3(ireg,imod,isce) = Changes (DGeo3(ireg,ifile), DGeo3(ireg,imod+1))

         ! Del_DDur1(ireg,imod,isce) = (DDur1(ireg,ifile)-DDur1(ireg,imod+1))/DDur1(ireg,imod+1)*100.
         ! Del_DDur2(ireg,imod,isce) = (DDur2(ireg,ifile)-DDur2(ireg,imod+1))/DDur2(ireg,imod+1)*100.
         ! Del_DDur3(ireg,imod,isce) = (DDur3(ireg,ifile)-DDur3(ireg,imod+1))/DDur3(ireg,imod+1)*100.

         Del_DDur1(ireg,imod,isce) = Changes (DDur1(ireg,ifile), DDur1(ireg,imod+1))
         Del_DDur2(ireg,imod,isce) = Changes (DDur2(ireg,ifile), DDur2(ireg,imod+1))
         Del_DDur3(ireg,imod,isce) = Changes (DDur3(ireg,ifile), DDur3(ireg,imod+1))

         ! Del_DInt1(ireg,imod,isce) = (DInt1(ireg,ifile)-DInt1(ireg,imod+1))/DInt1(ireg,imod+1)*100.
         ! Del_DInt2(ireg,imod,isce) = (DInt2(ireg,ifile)-DInt2(ireg,imod+1))/DInt2(ireg,imod+1)*100.
         ! Del_DInt3(ireg,imod,isce) = (DInt3(ireg,ifile)-DInt3(ireg,imod+1))/DInt3(ireg,imod+1)*100.

         Del_DInt1(ireg,imod,isce) = Changes (DInt1(ireg,ifile), DInt1(ireg,imod+1))
         Del_DInt2(ireg,imod,isce) = Changes (DInt2(ireg,ifile), DInt2(ireg,imod+1))
         Del_DInt3(ireg,imod,isce) = Changes (DInt3(ireg,ifile), DInt3(ireg,imod+1))

         ! Del_DSev1(ireg,imod,isce) = (DSev1(ireg,ifile)-DSev1(ireg,imod+1))/DSev1(ireg,imod+1)*100.
         ! Del_DSev2(ireg,imod,isce) = (DSev2(ireg,ifile)-DSev2(ireg,imod+1))/DSev2(ireg,imod+1)*100.
         ! Del_DSev3(ireg,imod,isce) = (DSev3(ireg,ifile)-DSev3(ireg,imod+1))/DSev3(ireg,imod+1)*100.

         Del_DSev1(ireg,imod,isce) = Changes (DSev1(ireg,ifile), DSev1(ireg,imod+1))
         Del_DSev2(ireg,imod,isce) = Changes (DSev2(ireg,ifile), DSev2(ireg,imod+1))
         Del_DSev3(ireg,imod,isce) = Changes (DSev3(ireg,ifile), DSev3(ireg,imod+1))

         ! Del_DSpl1(ireg,imod,isce) = (DSpl1(ireg,ifile)-DSpl1(ireg,imod+1))/DSpl1(ireg,imod+1)*100.
         ! Del_DSpl2(ireg,imod,isce) = (DSpl2(ireg,ifile)-DSpl2(ireg,imod+1))/DSpl2(ireg,imod+1)*100.
         ! Del_DSpl3(ireg,imod,isce) = (DSpl3(ireg,ifile)-DSpl3(ireg,imod+1))/DSpl3(ireg,imod+1)*100.

         Del_DSpl1(ireg,imod,isce) = Changes (DSpl1(ireg,ifile), DSpl1(ireg,imod+1)) 
         Del_DSpl2(ireg,imod,isce) = Changes (DSpl2(ireg,ifile), DSpl2(ireg,imod+1))
         Del_DSpl3(ireg,imod,isce) = Changes (DSpl3(ireg,ifile), DSpl3(ireg,imod+1))
      Enddo
   Enddo
Enddo

! Differences of Drought Characteristics between BL and OBS: May-03, 2018

Do ilon=1, Nlon
   Do ilat=1, Nlat
      if (Reg(ilon,ilat)>=1.AND.Reg(ilon,ilat)<=Nreg) then
         Do imod=1, Nmodels+1
            ifile=imod+1         ! ifile: Baseline file; 1: OBS 
            Del_BL_Num_of_Mon_D1(ilon,ilat,imod) = Changes (Nmon_D1(ilon,ilat,ifile), Nmon_D1(ilon,ilat,1))
            Del_BL_Num_of_Mon_D2(ilon,ilat,imod) = Changes (Nmon_D2(ilon,ilat,ifile), Nmon_D2(ilon,ilat,1))
            Del_BL_Num_of_Mon_D3(ilon,ilat,imod) = Changes (Nmon_D3(ilon,ilat,ifile), Nmon_D3(ilon,ilat,1))

            Del_BL_N_Spel_Grd_D1(ilon,ilat,imod) = Changes (NSpell_D1(ilon,ilat,ifile), NSpell_D1(ilon,ilat,1))
            Del_BL_N_Spel_Grd_D2(ilon,ilat,imod) = Changes (NSpell_D2(ilon,ilat,ifile), NSpell_D2(ilon,ilat,1))
            Del_BL_N_Spel_Grd_D3(ilon,ilat,imod) = Changes (NSpell_D3(ilon,ilat,ifile), NSpell_D3(ilon,ilat,1))

            Del_BL_Droug_Seve_D1(ilon,ilat,imod) = Changes (Sum_DS1(ilon,ilat,ifile), Sum_DS1(ilon,ilat,1))
            Del_BL_Droug_Seve_D2(ilon,ilat,imod) = Changes (Sum_DS2(ilon,ilat,ifile), Sum_DS2(ilon,ilat,1))
            Del_BL_Droug_Seve_D3(ilon,ilat,imod) = Changes (Sum_DS3(ilon,ilat,ifile), Sum_DS3(ilon,ilat,1))

            Del_BL_Mean_Durat_D1(ilon,ilat,imod) = Changes (Mean_Duration_D1(ilon,ilat,ifile), Mean_Duration_D1(ilon,ilat,1))
            Del_BL_Mean_Durat_D2(ilon,ilat,imod) = Changes (Mean_Duration_D2(ilon,ilat,ifile), Mean_Duration_D2(ilon,ilat,1))    
            Del_BL_Mean_Durat_D3(ilon,ilat,imod) = Changes (Mean_Duration_D3(ilon,ilat,ifile), Mean_Duration_D3(ilon,ilat,1))    
   
            Del_BL_Inter_Arri_D1(ilon,ilat,imod) = Changes (Inter_Arr_Time_D1(ilon,ilat,ifile), Inter_Arr_Time_D1(ilon,ilat,1))
            Del_BL_Inter_Arri_D2(ilon,ilat,imod) = Changes (Inter_Arr_Time_D2(ilon,ilat,ifile), Inter_Arr_Time_D2(ilon,ilat,1))
            Del_BL_Inter_Arri_D3(ilon,ilat,imod) = Changes (Inter_Arr_Time_D3(ilon,ilat,ifile), Inter_Arr_Time_D3(ilon,ilat,1))
         Enddo
      endif
   Enddo
Enddo
!
Do ireg=1, Nreg+1       
   Do imod=1, Nmodels+1
      ifile=imod+1         ! ifile: Baseline file; 1: OBS 
      Del_BL_Reg_Dr_Mon_D1(ireg,imod) = Changes (DMon1(ireg,ifile), DMon1(ireg,1))
      Del_BL_Reg_Dr_Mon_D2(ireg,imod) = Changes (DMon2(ireg,ifile), DMon2(ireg,1))
      Del_BL_Reg_Dr_Mon_D3(ireg,imod) = Changes (DMon3(ireg,ifile), DMon3(ireg,1))

      Del_BL_Reg_N_Spel_D1(ireg,imod) = Changes (DSpl1(ireg,ifile), DSpl1(ireg,1))
      Del_BL_Reg_N_Spel_D2(ireg,imod) = Changes (DSpl2(ireg,ifile), DSpl2(ireg,1))
      Del_BL_Reg_N_Spel_D3(ireg,imod) = Changes (DSpl3(ireg,ifile), DSpl3(ireg,1))

      Del_BL_Reg_Dr_Sev_D1(ireg,imod) = Changes (DSev1(ireg,ifile), DSev1(ireg,1))
      Del_BL_Reg_Dr_Sev_D2(ireg,imod) = Changes (DSev2(ireg,ifile), DSev2(ireg,1))
      Del_BL_Reg_Dr_Sev_D3(ireg,imod) = Changes (DSev3(ireg,ifile), DSev3(ireg,1))

      Del_BL_Reg_Durati_D1(ireg,imod) = Changes (DSev1(ireg,ifile), DSev1(ireg,1))
      Del_BL_Reg_Durati_D2(ireg,imod) = Changes (DSev2(ireg,ifile), DSev2(ireg,1))
      Del_BL_Reg_Durati_D3(ireg,imod) = Changes (DSev3(ireg,ifile), DSev3(ireg,1))

      Del_BL_Reg_Int_Ar_D1(ireg,imod) = Changes (DInt1(ireg,ifile), DInt1(ireg,1))
      Del_BL_Reg_Int_Ar_D2(ireg,imod) = Changes (DInt2(ireg,ifile), DInt2(ireg,1))
      Del_BL_Reg_Int_Ar_D3(ireg,imod) = Changes (DInt3(ireg,ifile), DInt3(ireg,1))

      Del_BL_Reg_Geo_Ex_D1(ireg,imod) = Changes (DGeo1(ireg,ifile), DGeo1(ireg,1))
      Del_BL_Reg_Geo_Ex_D2(ireg,imod) = Changes (DGeo2(ireg,ifile), DGeo2(ireg,1))
      Del_BL_Reg_Geo_Ex_D3(ireg,imod) = Changes (DGeo3(ireg,ifile), DGeo3(ireg,1))
   Enddo
Enddo

! Calculating the Agreement:

Call Agreement(DelNmon_D1, Agree_Nmon1); Call Agreement(DelNmon_D2, Agree_Nmon2); Call Agreement(DelNmon_D3, Agree_Nmon3)
Call Agreement(DelSum_DS1, Agree_DS1);   Call Agreement(DelSum_DS2, Agree_DS2);   Call Agreement(DelSum_DS3, Agree_DS3)
Call Agreement(Del_Mean_Duration_D1, Agree_Dur1);Call Agreement(Del_Mean_Duration_D2, Agree_Dur2);Call Agreement(Del_Mean_Duration_D3, Agree_Dur3)
Call Agreement(Del_Inter_Arr_Time_D1, Agree_Int1);Call Agreement(Del_Inter_Arr_Time_D2, Agree_Int2);Call Agreement(Del_Inter_Arr_Time_D3, Agree_Int3)

End Subroutine Calc_Drought_Characteristics
!------------------------------------------
!===========================================================================

Real Function Changes (Future, BL)
Real :: Future, BL, tmp
  if (Future /= Miss .AND. BL /= Miss .AND. BL /= 0) then
     tmp = (ABS(Future) - ABS(BL))/ABS(BL)*100.0
  else
     tmp = Miss
  endif
  Changes = tmp
End Function Changes
!------------------------------------------
Subroutine Ensemble_1 (X)
Integer :: ilon,ilat,ifile,imod,isce,it,ireg, imon,imon1, if1, if2
Real :: X(Nlon,Nlat,Nfile), Y(Nmodels)
! Ensamble of the Drought Characteristics:
Do ilon=1, Nlon
   Do ilat=1, Nlat
      Do isce=1, NSCs+1   ! 4 Scenarios + 1 BL
         if1=(isce-1)*(Nmodels+1) + 2; if2=if1+(Nmodels-1); Y(1:Nmodels) = X(ilon,ilat,if1:if2); X(ilon,ilat,if2+1)=Time_Ave(Y, Nmodels)
      Enddo
   Enddo
Enddo
End Subroutine Ensemble_1
!------------------------------------------
Subroutine Ensemble_2 (X)
Integer :: ilon,ilat,ifile,imod,isce,it,ireg, imon,imon1, if1, if2
Real :: X(Nreg+1,Nmon,Nfile), Y(Nmodels) ! Sum_Geo_D1(ireg,imon,ifile)
! Ensamble of the Drought Characteristics:
Do ireg=1, Nreg+1
   Do imon=1, Nmon
      Do isce=1, NSCs+1   ! 4 Scenarios + 1 BL
         if1=(isce-1)*(Nmodels+1) + 2; if2=if1+(Nmodels-1); Y(1:Nmodels) = X(ireg,imon,if1:if2); X(ireg,imon,if2+1)=Time_Ave(Y, Nmodels)
      Enddo
   Enddo
Enddo
End Subroutine Ensemble_2
!------------------------------------------
Subroutine Agreement (X, Y)
Integer :: ilon, ilat, isce, imod
Real :: X(Nlon,Nlat,Nmodels+1,NSCs), Y(Nlon,Nlat,NSCs), ND1, ASign
Y(:,:,:)=Miss
Do ilon=1, Nlon
   Do ilat=1, Nlat
   if (Reg(ilon,ilat)>0) then
      Do isce=1, NSCs
         ND1=0; ASign=X(ilon,ilat,Nmodels+1,isce)/ABS(X(ilon,ilat,Nmodels+1,isce))
         Do imod=1, Nmodels
            if (X(ilon,ilat,imod,isce)*X(ilon,ilat,Nmodels+1,isce)>0) ND1=ND1+1 
         Enddo
         Y(ilon,ilat,isce)=ND1/Real(Nmodels)*ASign*100.
      Enddo
   endif
   Enddo
Enddo 
End Subroutine Agreement
!------------------------------------------
Subroutine Regional_Mean (X, Y)
Integer :: ifile, ilon, ilat, ireg, Nc
Real, Dimension (Nlon,Nlat,Nfile) :: X
Real, dimension (Nreg+1,Nfile)    :: Y
Real :: tmp
Do ifile=1, Nfile
   Do ireg=1, Nreg+1
      if (ireg>=1.AND.ireg<=Nreg) then
         tmp=0; Nc=0
         Do ilon=1, Nlon
            Do ilat=1, Nlat
               if (Int(Reg(ilon,ilat))==ireg) then
                  if (X(ilon,ilat,ifile)/=Miss.AND.X(ilon,ilat,ifile)/=0.) then; Nc=Nc+1; tmp=tmp+X(ilon,ilat,ifile); endif
               endif
            Enddo
         Enddo
         if (Nc>1) then; Y(ireg,ifile)=tmp/Real(Nc); else; Y(ireg,ifile)=Miss; endif
      else if (ireg==Nreg+1) then
         tmp=0; Nc=0
         Do ilon=1, Nlon
            Do ilat=1, Nlat
               if (Int(Reg(ilon,ilat))>=1.AND.Int(Reg(ilon,ilat))<=Nreg) then
                  if (X(ilon,ilat,ifile)/=Miss.AND.X(ilon,ilat,ifile)/=0.) then; Nc=Nc+1; tmp=tmp+X(ilon,ilat,ifile); endif
               endif
            Enddo
         Enddo
         if (Nc>1) then; Y(ireg,ifile)=tmp/Real(Nc); else; Y(ireg,ifile)=Miss; endif
      endif
   Enddo
Enddo
End Subroutine Regional_Mean
!------------------------------------------