Subroutine Calc_Drought_Characteristics

Integer :: ist,ifile,imod,isce,it,ireg, imon,imon1, if1, ienso
Real :: ND1, ND2, ND3, Ngrd, NEL, N1, N2, N3
Real :: X(Nmon), XX(Nst), YY(Nreg+1)

! Frequency (NDmon/NEnsoMon) with respect to ENSO phases
  Do ist=1, Nst
     Do ienso=1, NEnso
        ND1=0; ND2=0; ND3=0; NEL=0
        Do it=1, Ntime
           if (ENSO(it)==ENSO_Phase(ienso)) then
              NEL=NEL+1
              if (D1(ist,it)==1) ND1=ND1+1; if (D2(ist,it)==1) ND2=ND2+1; if (D3(ist,it)==1) ND3=ND3+1
           endif
        Enddo
        Fre_ENSO_D1(ist,ienso)=ND1/NEL*100.; Fre_ENSO_D2(ist,ienso)=ND2/NEL*100.; Fre_ENSO_D3(ist,ienso)=ND3/NEL*100.
     Enddo
  Enddo
  ! Regional mean:
  Do ienso=1, NEnso
     XX(1:Nst)=Fre_ENSO_D1(1:Nst,ienso); Call Regional_Mean (XX, YY); Fre_ENSO_Reg_D1(1:Nreg+1,ienso)=YY(1:Nreg+1)
     XX(1:Nst)=Fre_ENSO_D2(1:Nst,ienso); Call Regional_Mean (XX, YY); Fre_ENSO_Reg_D2(1:Nreg+1,ienso)=YY(1:Nreg+1)
     XX(1:Nst)=Fre_ENSO_D3(1:Nst,ienso); Call Regional_Mean (XX, YY); Fre_ENSO_Reg_D3(1:Nreg+1,ienso)=YY(1:Nreg+1)
  Enddo
! Severity with respect to ENSO phases = Sum(PDSI)/Sum(ENSO_mon)
  Do ist=1, Nst
     Do ienso=1, NEnso
        ND1=0; ND2=0; ND3=0; NEL=0; N1=0; N2=0; N3=0
        Do it=1, Ntime
           if (ENSO(it)==ENSO_Phase(ienso)) then
              NEL=NEL+1
!              if (D1(ist,it)==1) ND1=ND1+PDSI(ist,it); if (D2(ist,it)==1) ND2=ND2+PDSI(ist,it); if (D3(ist,it)==1) ND3=ND3+PDSI(ist,it)
              if (D1(ist,it)==1) then; N1=N1+1; ND1=ND1+PDSI(ist,it); endif
              if (D2(ist,it)==1) then; N2=N2+1; ND2=ND2+PDSI(ist,it); endif
              if (D3(ist,it)==1) then; N3=N3+1; ND3=ND3+PDSI(ist,it); endif
           endif
        Enddo
!        Sev_ENSO_D1(ist,ienso)=ND1/NEL; Sev_ENSO_D2(ist,ienso)=ND2/NEL; Sev_ENSO_D3(ist,ienso)=ND3/NEL
!        Sev_ENSO_D1(ist,ienso)=ND1/N1; Sev_ENSO_D2(ist,ienso)=ND2/N2; Sev_ENSO_D3(ist,ienso)=ND3/N3
        if (N1>0) Sev_ENSO_D1(ist,ienso)=ND1/N1
        if (N2>0) Sev_ENSO_D2(ist,ienso)=ND2/N2
        if (N3>0) Sev_ENSO_D3(ist,ienso)=ND3/N3
     Enddo
  Enddo
  ! Regional mean:
  Do ienso=1, NEnso
     XX(1:Nst)=Sev_ENSO_D1(1:Nst,ienso); Call Regional_Mean (XX, YY); Sev_ENSO_Reg_D1(1:Nreg+1,ienso)=YY(1:Nreg+1)
     XX(1:Nst)=Sev_ENSO_D2(1:Nst,ienso); Call Regional_Mean (XX, YY); Sev_ENSO_Reg_D2(1:Nreg+1,ienso)=YY(1:Nreg+1)
     XX(1:Nst)=Sev_ENSO_D3(1:Nst,ienso); Call Regional_Mean (XX, YY); Sev_ENSO_Reg_D3(1:Nreg+1,ienso)=YY(1:Nreg+1)
  Enddo

! Geo Extent ENSO:
  Do ireg=1, Nreg+1
     Do ienso=1, NEnso
        ND1=0; ND2=0; ND3=0; NEL=0; N1=0; N2=0; N3=0
        Do it=1, Ntime
           if (ENSO(it)==ENSO_Phase(ienso)) then  ! ENSO phase
              NEL=NEL+1
              if (Geo_D1(ireg,it)/=Miss) then; N1=N1+1; ND1=ND1+Geo_D1(ireg,it); endif
              if (Geo_D2(ireg,it)/=Miss) then; N2=N2+1; ND2=ND2+Geo_D2(ireg,it); endif
              if (Geo_D3(ireg,it)/=Miss) then; N3=N3+1; ND3=ND3+Geo_D3(ireg,it); endif
           endif
!           Geo_ENSO_Reg_D1(ireg,ienso)=ND1/NEL; Geo_ENSO_Reg_D2(ireg,ienso)=ND2/NEL; Geo_ENSO_Reg_D3(ireg,ienso)=ND3/NEL;
           Geo_ENSO_Reg_D1(ireg,ienso)=ND1/N1; Geo_ENSO_Reg_D2(ireg,ienso)=ND2/N2; Geo_ENSO_Reg_D3(ireg,ienso)=ND3/N3;
        Enddo 
     Enddo
  Enddo 

! Regional mean of T2m & Prec:
   Do it=1, Ntime
      XX(1:Nst)=T2m (1:Nst,it); Call Regional_Mean (XX, YY); T2m_Reg(1:Nreg+1,it)=YY(1:Nreg+1)
      XX(1:Nst)=Prec(1:Nst,it); Call Regional_Mean (XX, YY); Pre_Reg(1:Nreg+1,it)=YY(1:Nreg+1)
   Enddo

   Do ireg=1, Nreg+1
      Do imon=1, Nmon
         ND1=0; ND2=0; 
         Do it=1, Ntime
            imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon
            if (imon1==imon) then
               ND1=ND1+T2m_Reg(ireg,it); ND2=ND2+Pre_Reg(ireg,it)
            endif
         Enddo
         T2m_Reg_Mon(ireg,imon)=ND1/Real(NYear); Pre_Reg_Mon(ireg,imon)=ND2/Real(NYear)
      Enddo
   Enddo

! Sum over Ntime of the Number of Drought Month D1,2,3
   Do ist=1, Nst
         ND1=0; ND2=0; ND3=0
         Do it=1, Ntime
            if (D1(ist,it)==1.) ND1=ND1+1
            if (D2(ist,it)==1.) ND2=ND2+1
            if (D3(ist,it)==1.) ND3=ND3+1
         Enddo
         Nmon_D1(ist)=ND1/Ntime; Nmon_D2(ist)=ND2/Ntime; Nmon_D3(ist)=ND3/Ntime   ! Drought Frequency = Num of Dr mon / Ntime
   Enddo

! Monthly Num of Drou Mon at each station (Nst, Nmon) = Freq 
   Do ist=1, Nst
      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
               if (D1(ist,it)==1.) ND1=ND1+D1(ist,it); 
               if (D2(ist,it)==1.) ND2=ND2+D2(ist,it); 
               if (D3(ist,it)==1.) ND3=ND3+D3(ist,it);
            endif
         Enddo
         Mon_NDM_D1(ist,imon)=ND1/Real(NYear); Mon_NDM_D2(ist,imon)=ND2/Real(NYear); Mon_NDM_D3(ist,imon)=ND3/Real(NYear);
         ! Average over all year
      Enddo
   Enddo

! Monthly Severity at each station (Nst, Nmon) 
   Do ist=1, Nst
      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
               if (D1(ist,it)==1.) ND1=ND1+PDSI(ist,it)
               if (D2(ist,it)==1.) ND2=ND2+PDSI(ist,it)
               if (D3(ist,it)==1.) ND3=ND3+PDSI(ist,it)
            endif
         Enddo
         Mon_Sev_D1(ist,imon)=ND1; Mon_Sev_D2(ist,imon)=ND2; Mon_Sev_D3(ist,imon)=ND3;
         ! Mon_Sev_D1(ist,imon)=ND1/Real(NYear); Mon_Sev_D2(ist,imon)=ND2/Real(NYear); Mon_Sev_D3(ist,imon)=ND3/Real(NYear);
         ! Average over all year
      Enddo
   Enddo

! Sum over Ntime of Drought Severity
   Do ist=1, Nst
         ND1=0; ND2=0; ND3=0
         Do it=1, Ntime
            if (D1(ist,it)==1.) ND1=ND1+PDSI(ist,it)
            if (D2(ist,it)==1.) ND2=ND2+PDSI(ist,it)
            if (D3(ist,it)==1.) ND3=ND3+PDSI(ist,it)
         Enddo
         Sum_DS1(ist)=ND1; Sum_DS2(ist)=ND2; Sum_DS3(ist)=ND3
   Enddo

! Monthly mean of Geo_Ext = Sum over Ntime of Geo_D1,2,3 divide to Ntime
   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); ND2=ND2+Geo_D2(ireg,it); ND3=ND3+Geo_D3(ireg,it);
            endif
         Enddo
         Sum_Geo_D1(ireg,imon)=ND1/Real(NYear); Sum_Geo_D2(ireg,imon)=ND2/Real(NYear); Sum_Geo_D3(ireg,imon)=ND3/Real(NYear)
      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)

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

Do imon=1, Nmon
!  XX(Nst), YY(Nreg+1)
   XX(1:Nst)=Mon_NDM_D1(1:Nst,imon); Call Regional_Mean (XX, YY); Sum_NDM_D1(1:Nreg+1,imon)=YY(1:Nreg+1)
   XX(1:Nst)=Mon_NDM_D2(1:Nst,imon); Call Regional_Mean (XX, YY); Sum_NDM_D2(1:Nreg+1,imon)=YY(1:Nreg+1)
   XX(1:Nst)=Mon_NDM_D3(1:Nst,imon); Call Regional_Mean (XX, YY); Sum_NDM_D3(1:Nreg+1,imon)=YY(1:Nreg+1)

   XX(1:Nst)=Mon_Sev_D1(1:Nst,imon); Call Regional_Mean (XX, YY); Sum_Sev_D1(1:Nreg+1,imon)=YY(1:Nreg+1)
   XX(1:Nst)=Mon_Sev_D2(1:Nst,imon); Call Regional_Mean (XX, YY); Sum_Sev_D2(1:Nreg+1,imon)=YY(1:Nreg+1)
   XX(1:Nst)=Mon_Sev_D3(1:Nst,imon); Call Regional_Mean (XX, YY); Sum_Sev_D3(1:Nreg+1,imon)=YY(1:Nreg+1)   
Enddo

End Subroutine Calc_Drought_Characteristics

!------------------------------------------
Subroutine Regional_Mean (X, Y)
Integer :: ifile, ist, ireg, Nc
Real, Dimension (Nst) :: X
Real, dimension (Nreg+1)    :: Y
Real :: tmp
   Do ireg=1, Nreg+1
      if (ireg>=1.AND.ireg<=Nreg) then
         tmp=0; Nc=0
         Do ist=1, Nst
               if (Int(Reg(ist))==ireg) then
                  if (X(ist)/=Miss.AND.X(ist)/=0.) then; Nc=Nc+1; tmp=tmp+X(ist); endif
               endif
         Enddo
         if (Nc>1) then; Y(ireg)=tmp/Real(Nc); else; Y(ireg)=Miss; endif
      else if (ireg==Nreg+1) then
         tmp=0; Nc=0
         Do ist=1, Nst
               if (Int(Reg(ist))>=1.AND.Int(Reg(ist))<=Nreg) then
                  if (X(ist)/=Miss.AND.X(ist)/=0.) then; Nc=Nc+1; tmp=tmp+X(ist); endif
               endif
         Enddo
         if (Nc>1) then; Y(ireg)=tmp/Real(Nc); else; Y(ireg)=Miss; endif
      endif
   Enddo
End Subroutine Regional_Mean
!------------------------------------------