      subroutine PDINEW(iendyr,iendmo,IDAY,IMULT,DAYS)
c Subroutine to compute Palmer Drought Severity Index (PDSI), Palmer 
c Hydrologic Drought Index (PHDI), real-time Palmer Drought Index (WPLM),
c and Z index
c Original program written by Ned Guttman and modified by Richard Heim.
c Extensively modified to run as a subroutine on shared drive,
c March 2016, Trevor Wallis
c Modified to output a current PHDI file for drought ameliorarion
c project, Trevor wallis April 2017
c23456789012345678901234567890123456789012345678901234567890123456789012
c
       COMMON /COMPATHS/ PATHCPC,PATHDRD,PATHDRDW,cMeanpath,PATHWDATA,
     :                   cOdata,cIndata
       CHARACTER PATHCPC*42,PATHDRD*34,PATHDRDW*40,cMeanpath*36,
     :           PATHWDATA*41,cIndata*38,cOdata*38,name*20,
     :           cName(10:37)*20,cI*2,cEl(4)*4,cSt*2
       integer ISTAT(12,4),DAYS(12),iOut(10:37),indexj(40),indexm(40),
     :         lOut(10:37)
       real PHI(12),prec(200,14),temp(200,14),prec1(200,14),
     :      temp1(200,14),pdat(200,12),spdat(200,12),pedat(200,12),
     :      pldat(200,12),prdat(200,12),rdat(200,12),tldat(200,12),
     :      etdat(200,12),rodat(200,12),tdat(200,12),sssdat(200,12),
     :      ssudat(200,12),tsum(12),psum(12),spsum(12),pesum(12),
     :      plsum(12),prsum(12),rsum(12),tlsum(12),etsum(12),rosum(12),
     *      pdsi(200,13),phdi(200,13),wplm(200,13),alp(12),bet(12),
     *      gam(12),del(12),trat(12),sabsd(12),akhat(12),ak(12),
     *      cp(200,13),z(200,13),ppr(200,12),px1(200,12),
     *      px2(200,12),px3(200,12),eff(200,12),
     *      sx(96),sx1(40),sx2(40),sx3(40),x(200,12),phatsav(12)

      DATA PHI  /-.3865982,-.2316132,-.0378180, .1715539, .3458803, 
     *            .4308320, .3916645, .2452467, .0535511,-.15583436,
     *           -.3340551,-.4310691 /  
      data cEl /'pdsi','phdi','zndx','pmdi' /

C---GET STATUS INFO FOR PUERTO RICO & OTHER TERRITORIES
       OPEN (4,FILE=PATHDRDW//'prolog.status',STATUS='OLD')
       IYRPR=1998
       LMNPR=6
401    READ (4,404,END=499) IY,((ISTAT(M,I),I=1,4),M=1,12)
404    FORMAT (1X,I4,12(1X,4I1))
       DO 410 M=1,12
       IF (ISTAT(M,2).EQ.0 .AND. ISTAT(M,3).EQ.0 .AND. 
     *     ISTAT(M,4).EQ.0) GO TO 499
       IYRPR=IY
       LMNPR=M
410    CONTINUE
       GO TO 401
499    CLOSE (4)
c Open parameter file
      open (7, file=cIndata//'pdinew.parameter',status='old')
c input parameters from pdinew.parameter file (unit 7)
      read (7,'(I4)')   nbcaly
       jbcaly=nbcaly
      read (7,'(I4)')   necaly
      read (7,'(I4)')   nbegyr
       jbegyr=nbegyr
      read (7,'(I4)')   ioutbyr
       joutbyr=ioutbyr
      read (7,'(I4)')   iouteyr
      iouteyr=iendyr
      read (7,'(F6.2)') amsng
      read (7,'(I2)')   inunit
      close (7)
c Read and open output file names
      open(7,file=PATHDRDW//'pdinew.outputfiles',status='old')
      do 101,i=1,24
       read(7,7,end=1007) name
 7     format(a20)
       if ((i .ge. 1) .and.(i .le. 6)) then
        iOut(i+9) = inunit + (i - 1) + 10
        cName(i+9) = name
        open(iOut(i+9),file=cOdata//'palmer/'//cName(i+9),
     :         status='unknown')
       else if (i .ge. 7) then
        iOut(i+13) = inunit + (i - 1) + 13
        cName(i+13) = name
        open(iOut(i+13),file=cOdata//'palmer/'//cName(i+13),
     :         status='unknown')
       endif
 101  continue
 1007 close (7)
c Provide readout as to which output files are being processed
      if (cName(10)(13:17) .eq. 'div  ') then
       write(6,*) 'Processing climate division Palmer indices'
       open(425,file=cOdata//'amelioration/wmp-phdi.dat',
     :        status='unknown')
       write(425,425) 'Weekly-monthly Palmer for month to ',
     :        iendyr,iendmo,IDAY
 425   format(a35,i4.4,2i2.2)
       do 599,i=1,4
        write(cI,'(i2.2)') i + 4
        open(i+74,file=cOdata//'palmer/'//cEl(i)//cI//'dv.',
     :              status='old')
 599   continue
      else if (cName(10)(13:17) .eq. 'state') then
       write(6,*) 'Processing state and region Palmer indices'
       do 598,i=1,4
        write(cI,'(i2.2)') i + 4
        open(i+74,file=cOdata//'palmer/'//cEl(i)//cI//'st.',
     :              status='old')
 598   continue
      else if (cName(10)(13:17) .eq. 'river') then
       write(6,*) 'Processing river basin Palmer indices'
       do 597,i=1,4
        write(cI,'(i2.2)') i + 4
        open(i+74,file=cOdata//'palmer/'//cEl(i)//cI//'st.',
     :              status='old',access='append')
 597   continue
      endif
c---------------------------------------------------------------
c     read in soil constants
      open (114,file=PATHDRDW//'pdinew.soilconst',status='old')
 15   read(114,1002,END=16)IST,IDV,wcbot,b,h,tla
 1002 format(i3,i2,4f8.4) 
C---STATES ALASKA TO THE PACIFIC ISLANDS DO NOT HAVE DATA PRIOR
C---TO 1931.  SO, NEED TO SET START YEARS TO 1931 FOR THOSE
C---STATE CODES.
       IF (IST.GT.48 .AND. IST.LT.100) THEN
         IF (nbcaly.lt.1931) nbcaly=1931
         if (nbegyr.lt.1931) nbegyr=1931
         if (ioutbyr.lt.1931) ioutbyr=1931
         nendyr=IYRPR
         nendmo=LMNPR
       ELSE
         nbcaly=jbcaly
         nbegyr=jbegyr
         ioutbyr=joutbyr
         nendyr=iendyr
         nendmo=iendmo
       ENDIF
       IF (IST.EQ.1 .AND. IDV.EQ.1) THEN
         WRITE (9,*) 'processing divisions' 
       ELSE IF (IST.EQ.1 .AND. IDV.EQ.0) THEN
         WRITE (9,*) 'processing states & T.K. regions'
       ELSE IF (IST.EQ.111 .AND. IDV.EQ.0) THEN
         WRITE (9,*) 'processing river basins and special regions'
       ENDIF
c
 17   WCTOP = 1.0
      SS  = WCTOP
      SU = WCBOT
      WCTOT = WCBOT + WCTOP
      DO 20 M = 1,12
        TSUM  (M) = 0.
        PSUM  (M) = 0.
        SPSUM (M) = 0.
        PESUM (M) = 0.
        PLSUM (M) = 0.
        PRSUM (M) = 0.
        RSUM  (M) = 0.
        TLSUM (M) = 0.
        ETSUM (M) = 0.
        ROSUM (M) = 0.
        do j=1,nendyr-nbegyr+1
          pdsi(j,m)=amsng
          phdi(j,m)=amsng
          z(j,m)=amsng
          wplm(j,m)=amsng
          cp(j,m)=amsng
          pdat(j,m)=amsng
          spdat(j,m)=amsng
          pedat(j,m)=amsng
          pldat(j,m)=amsng
          prdat(j,m)=amsng
          rdat(j,m)=amsng
          tldat(j,m)=amsng
          etdat(j,m)=amsng
          rodat(j,m)=amsng
          tdat(j,m)=amsng
          sssdat(j,m)=amsng
          ssudat(j,m)=amsng
          eff(j,m)=0.0
          PPR(j,m)=0.0
          PX1(j,m)=0.0
          PX2(j,m)=0.0
          PX3(j,m)=0.0
        enddo
 20   continue
c Get division, state or region temperature and precipitation
      call GETPOR(IST,IDV,1,prec1)
      call GETPOR(IST,IDV,2,temp1)
      in=1
      do 221 k=1,200
       iy=INT(prec1(k,14))
       if (iy.lt.nbegyr) go to 221
        do m=1,14
         prec(in,m)=prec1(k,m)
         temp(in,m)=temp1(k,m)
        enddo
       in=in+1
  221 continue
c loop on years (j) and months (m)
      do 70  j=1,nendyr-nbegyr+1
       iy=INT(prec(j,14))
c IS IT A LEAP YEAR?
        LYEAR = (IY/4*4)/IY
        do 60 m=1,12
          if ((j.eq.nendyr-nbegyr+1).and.(m.gt.nendmo)) goto 60
          t=temp(j,m)
          p=prec(j,m)
C START THE WATER BALANCE CALCULATIONS
      SP = SS + SU
      PR = WCBOT + WCTOP - SP
C CALCULATE PE (POTENTIAL EVAPOTRANSPIRATION)   
      IF (T.LE.32.) THEN
          PE   = 0.0
      ELSE  
          DUM  = PHI(M)*TLA 
          DK   = ATAN(SQRT(1. - DUM*DUM)/DUM)   
          IF (DK.LT.0.) DK = 3.141593 + DK  
          DK   = (DK + .0157)/1.57  
          IF (T.GE.80.) THEN
              PE = (SIN(T/57.3 - .166) - .76)*DK
          ELSE  
              DUM=log(t-32)
              PE = (EXP(-3.863233 + B*1.715598 -B*log(h) + B*DUM))*DK 
          ENDIF 
      ENDIF 
C CONVERT DAILY TO MONTHLY
      IF (M.EQ.2.AND.LYEAR.EQ.1) THEN   
          PE = PE * 29. 
      ELSE  
          PE = PE * DAYS(M) 
      ENDIF 
C PL  POTENTIAL LOSS
      IF (SS.GE.PE) THEN
          PL  = PE  
      ELSE  
          PL   = ((PE - SS) * SU) / (WCBOT + WCTOP) + SS   
          PL   = AMIN1(PL,SP)   
      ENDIF 
C CALCULATE RECHARGE, RUNOFF, RESIDUAL MOISTURE, LOSS TO BOTH   
C SURFACE AND UNDER LAYERS, DEPENDING ON STARTING MOISTURE  
C CONTENT AND VALUES OF PRECIPITATION AND EVAPORATION.  
      IF (P.GE.PE) THEN 
C PRECIP EXCEEDS POTENTIAL EVAPORATION
          ET = PE   
          TL = 0.0  
          IF ((P - PE).GT.(WCTOP - SS)) THEN   
C EXCESS PRECIP RECHARGES UNDER LAYER AS WELL AS UPPER   
              RS = WCTOP - SS  
              SSS = WCTOP  
              IF ((P - PE -RS).LT.(WCBOT - SU)) THEN
C BOTH LAYERS CAN TAKE THE ENTIRE EXCESS  
                  RU = P - PE - RS  
                  RO = 0.0  
              ELSE  
C SOME RUNOFF OCCURS 
                  RU = WCBOT - SU   
                  RO = P - PE - RS - RU 
              ENDIF 
              SSU = SU + RU 
              R   = RS + RU 
          ELSE  
C ONLY TOP LAYER RECHARGED   
              R   = P - PE  
              SSS = SS + P - PE 
              SSU = SU  
              RO  = 0.0 
          ENDIF 
      ELSE  
C EVAPORATION EXCEEDS PRECIPITATION   
          R  = 0.0  
          IF (SS.GE.(PE - P)) THEN  
C EVAP FROM SURFACE LAYER ONLY  
              SL  = PE - P  
              SSS = SS - SL 
              UL  = 0.0 
              SSU = SU  
          ELSE  
C EVAP FROM BOTH LAYERS 
              SL  = SS  
              SSS = 0.0 
c Error found in UL calculation. Rectified 4/28/2017, Trevor Wallis
              UL  = (PE - P - SL) * SU / (WCTOT)
*              UL  = (PE - P - SL) * SU / (WCTOT + 1.0)  
              UL  = AMIN1(UL,SU)
              SSU = SU - UL 
          ENDIF 
          TL  = SL + UL 
          RO  = 0.0 
          ET  = P  + SL + UL
      ENDIF 
C FOR CALIBRATION YEARS, SUM VALUES NEEDED TO CALCULATE THE 
C NORMAL CLIMATE COEFFICIENTS (ALPHA, BETA, ETC.)   
      IF(IY.GE.NBCALY.AND.IY.LE.NECALY) THEN
          TSUM (M) = TSUM (M) + T   
          PSUM (M) = PSUM (M) + P   
          SPSUM(M) = SPSUM(M) + SP  
          PESUM(M) = PESUM(M) + PE  
          PLSUM(M) = PLSUM(M) + PL  
          PRSUM(M) = PRSUM(M) + PR  
          RSUM (M) = RSUM (M) + R   
          TLSUM(M) = TLSUM(M) + TL  
          ETSUM(M) = ETSUM(M) + ET  
          ROSUM(M) = ROSUM(M) + RO  
      ENDIF 
      pdat(j,m)=p
      spdat(j,m)=sp
      pedat(j,m)=pe
      pldat(j,m)=pl
      prdat(j,m)=pr
      rdat(j,m)=r
      tldat(j,m)=tl
      etdat(j,m)=et
      rodat(j,m)=ro
      tdat(j,m)=t
      sssdat(j,m)=sss
      ssudat(j,m)=ssu
      
      SS       = SSS
      SU       = SSU

   60 CONTINUE  
C END OF MONTH LOOP 
c Output the intermediate variables
      if(iy.ge.ioutbyr.and.iy.le.iouteyr)then
       write(iout(20),43)IST,IDV,iy,(tdat(j,m),m=1,12)
       write(iout(21),43)IST,IDV,iy,(pdat(j,m),m=1,12)
       write(iout(22),43)IST,IDV,iy,(spdat(j,m),m=1,12)
       write(iout(23),43)IST,IDV,iy,(pedat(j,m),m=1,12)
       write(iout(24),43)IST,IDV,iy,(pldat(j,m),m=1,12)
       write(iout(25),43)IST,IDV,iy,(prdat(j,m),m=1,12)
       write(iout(26),43)IST,IDV,iy,(rdat(j,m),m=1,12)
       write(iout(27),43)IST,IDV,iy,(tldat(j,m),m=1,12)
       write(iout(28),43)IST,IDV,iy,(etdat(j,m),m=1,12)
       write(iout(29),43)IST,IDV,iy,(rodat(j,m),m=1,12)
       write(iout(30),43)IST,IDV,iy,(sssdat(j,m),m=1,12)
       write(iout(31),43)IST,IDV,iy,(ssudat(j,m),m=1,12)
      endif
43    format (I3.3,I2.2,I5,13F7.2)
   70 CONTINUE  
C END OF YEAR LOOP
C LOOP 90 IS BY MONTH: CAFEC COEFFICIENTS CALCULATIONS AND  
C 'T' RATIO FOR EACH MONTH  
      DO 90 M = 1, 12   
C ALPHA CALCULATION
      IF(PESUM(M).NE.0.) THEN   
          ALP(M) = ETSUM(M) / PESUM(M)  
      ELSE  
          IF (ETSUM(M).EQ.0.) THEN  
              ALP(M) = 1.0  
          ELSE  
              ALP(M) = 0.0  
          ENDIF 
      ENDIF 
C BETA CALCULATION
      IF (PRSUM(M).NE.0.) THEN  
          BET(M) = RSUM(M)/PRSUM(M) 
      ELSE  
          IF (RSUM(M).EQ.0.) THEN   
              BET(M) = 1.0  
          ELSE  
              BET(M) = 0.0  
          ENDIF 
      ENDIF 
C GAMMA CALCULATION 
      IF (SPSUM(M).NE.0.) THEN  
          GAM(M) = ROSUM(M) / SPSUM(M)  
      ELSE  
          IF (ROSUM(M).EQ.0.) THEN  
              GAM(M) = 1.   
          ELSE  
              GAM(M) = 0.0  
          ENDIF 
      ENDIF 
C DELTA CALCULATION 
      IF (PLSUM(M).NE.0.) THEN  
          DEL(M) = TLSUM(M) / PLSUM(M)  
      ELSE  
          DEL(M) = 0.0  
      ENDIF 
C CALCULATION OF Z-INDEX WEIGHTING FACTORS (VARIABLE AK)
C TRAT(M) IS THE 'T' RATIO OF AVERAGE MOISTURE DEMAND TO AVERAGE
C MOISTURE SUPPLY IN MONTH M.   
      TRAT(M) = (PESUM(M)+RSUM(M)+ROSUM(M)) / (PSUM(M)+TLSUM(M))
C 'SUM' VARIABLES NOW BECOME AVERAGES OVER THE CALIBRATION PERIOD   
      numyr=necaly-nbcaly+1
      TSUM (M) = TSUM(M)  / NUMYR   
      PSUM (M) = PSUM (M) / NUMYR   
      SPSUM(M) = SPSUM(M) / NUMYR   
      PESUM(M) = PESUM(M) / NUMYR   
      PLSUM(M) = PLSUM(M) / NUMYR   
      PRSUM(M) = PRSUM(M) / NUMYR   
      RSUM (M) = RSUM (M) / NUMYR   
      TLSUM(M) = TLSUM(M) / NUMYR   
      ETSUM(M) = ETSUM(M) / NUMYR   
      ROSUM(M) = ROSUM(M) / NUMYR   
   90 CONTINUE  
C END OF MONTHLY COEFFICIENT LOOP
      DO 100,M = 1, 12  
      SABSD(M)  = 0.0   
  100 continue          
C REREAD MONTHLY PARAMETERS FOR CALCULATION OF 
C THE 'K' MONTHLY WEIGHTING FACTORS USED IN Z-INDEX CALCULATION 
      do 102,j=nbcaly-nbegyr+1,necaly-nbegyr+1
       iy=j+nbegyr-1
*      do 102,j=nbegyr,nendyr
*       iy = j
       do 601,m=1,12
        PHAT = ALP(M)*PEdat(j,m) + BET(M)*PRdat(j,m) +
     *           GAM(M)*SPdat(j,m) - DEL(M)*PLdat(j,m)  
        phatsav(m)=PHAT
*        if ((j .ge. nbcaly-nbegyr+1) .and.
*     :      (j .le. necaly-nbegyr+1)) then
         D = Pdat(j,m) - PHAT   
         SABSD(M)= SABSD(M) + ABS(D) 
*        endif
 601   continue
c---save the PHAT values
       write (iout(37),43) IST,IDV,iy,(phatsav(m),m=1,12)
 102  continue
      SWTD = 0.0
      DO 130 MO = 1,12 
       DBAR = SABSD(MO) / NUMYR 
       AKHAT(MO) = 1.5 * log10((trat(mo)+2.8)/dbar) + .5
       SWTD = SWTD + DBAR*AKHAT(MO)  
  130 CONTINUE  
      DO 140 MO = 1,12  
       AK(MO) = 17.67 * AKHAT(MO) / SWTD 
C PRINT OUT COEFFICIENTS AND OTHER STUFF for calibration years
       write(iout(10),6060) IST,IDV,MO,ALP(MO),BET(MO),GAM(MO), 
     *   DEL(MO),AK(MO),PSUM(MO),SPSUM(MO),PRSUM(MO),PESUM(MO), 
     *   PLSUM(MO), ETSUM(MO),RSUM(MO),TLSUM(MO),ROSUM(MO)  
 140  continue
 6060 format(i3.3,i2.2,i3,4f9.4,f9.3,9f9.2)
C INITIALIZE VARIABLES
      V   = 0.0 
      PRO = 0.0 
      X1  = 0.0 
      X2  = 0.0 
      X3  = 0.0 
      K8  = 1
      k8max = 1
C REREAD data FOR CALCULATION OF Z-INDEX (MOISTURE ANOMALY) AND PDSI
C (VARIABLE X). FINAL OUTPUTS ARE VARIABLES PX3, X, AND Z
      do 231 j=1,nendyr-nbegyr+1
       iy=j+nbegyr-1
       do 230 m=1,12
        if ((j.eq.nendyr-nbegyr+1).and.(m.gt.nendmo)) then
         goto 2311
        endif
        indexj(k8)=j
        indexm(k8)=m
  160   ZE  = 0.0 
        UD  = 0.0 
        UW  = 0.0 
        CET = ALP(M) * PEdat(j,m)
        CR  = BET(M) * PRdat(j,m)
        CRO = GAM(M) * SPdat(j,m)
        CL  = DEL(M) * PLdat(j,m)
        CP(j,m)  = CET + CR + CRO - CL 
        CD  = Pdat(j,m)- CP(j,m)  
        Z(j,m)   = AK(M) * CD 
        IF (PRO.EQ.100. .OR. PRO.EQ.0.) THEN  
C NO ABATEMENT UNDERWAY WET OR DROUGHT WILL END IF -.5 =< X3 =< .5   
          IF (X3.GE.-.5 .AND. X3.LE..5) THEN
C END OF DROUGHT OR WET  
              PV  = 0.0 
              PPR(j,m) = 0.0 
              PX3(j,m) = 0.0 
C BUT CHECK FOR NEW WET OR DROUGHT START FIRST 
              GO TO 200 
          ELSE IF (X3.GT..5) THEN   
C WE ARE IN A WET SPELL 
               IF (Z(j,m).GE..15) THEN   
C THE WET SPELL INTENSIFIES 
                  GO TO 210 
              ELSE  
C THE WET STARTS TO ABATE (AND MAY END)  
                  GO TO 170 
              ENDIF 
          ELSE IF (X3.LT.-.5) THEN  
C WE ARE IN A DROUGHT 
               IF (Z(j,m).LE.-.15) THEN  
C THE DROUGHT INTENSIFIES 
                  GO TO 210 
              ELSE  
C THE DROUGHT STARTS TO ABATE (AND MAY END)  
                  GO TO 180 
              ENDIF 
          ENDIF 
      ELSE  
C ABATEMENT IS UNDERWAY   
          IF (X3.GT.0.) GO TO 170   
C WE ARE IN A WET SPELL 
          IF (X3.LE.0.) GO TO 180   
C WE ARE IN A DROUGHT   
      ENDIF 
C WET SPELL ABATEMENT IS POSSIBLE  
  170 UD = Z(j,m) - .15  
      PV = UD + AMIN1(V,0.) 
      IF (PV.GE.0) GO TO 210
C DURING A WET SPELL, PV => 0 IMPLIES PROB(END) HAS RETURNED TO 0
      ZE = -2.691 * X3 + 1.5
      GO TO 190 
C DROUGHT ABATEMENT IS POSSIBLE
  180 UW = Z(j,m) + .15  
      PV = UW + AMAX1(V,0.) 
      IF (PV.LE.0) GO TO 210
C DURING A DROUGHT, PV =< 0 IMPLIES PROB(END) HAS RETURNED TO 0
      ZE = -2.691 * X3 - 1.5
C PROB(END) = 100 * (V/Q)  WHERE:   
C    V = SUM OF MOISTURE EXCESS OR DEFICIT (UD OR UW)  
C           DURING CURRENT ABATEMENT PERIOD, AND
C    Q = TOTAL MOISTURE ANOMALY REQUIRED TO END THE
C           CURRENT DROUGHT OR WET SPELL  
  190 IF (PRO.EQ.100.) THEN 
C DROUGHT OR WET CONTINUES, CALCULATE PROB(END) - VARIABLE ZE 
          Q = ZE
      ELSE  
          Q = ZE + V
      ENDIF 
      PPR(j,m) = (PV / Q) * 100. 
      IF (PPR(j,m).GE.100.) THEN 
          PPR(j,m) = 100.
          PX3(j,m) = 0.  
      ELSE  
          PX3(j,m) = .897 * X3 + Z(j,m)/3.
      ENDIF 
C CONTINUE X1 AND X2 CALCULATIONS. IF EITHER INDICATES THE START OF A NEW
C WET OR DROUGHT, AND IF THE LAST WET OR DROUGHT HAS ENDED, USE X1 OR X2
C AS THE NEW X3.
  200 PX1(j,m) = .897 * X1 + Z(j,m)/3.
      PX1(j,m) = AMAX1(PX1(j,m),0.)   
      IF (PX1(j,m).GE.1.) THEN   
          IF (PX3(j,m).EQ.0.) THEN   
C IF NO EXISTING WET SPELL OR DROUGHT X1 BECOMES THE NEW X3 
              X(j,m)   = PX1(j,m) 
              PX3(j,m) = PX1(j,m) 
              PX1(j,m) = 0. 
              iass=1
              CALL ASSIGN (iass,k8,ppr,px1,px2,px3,x,pdsi,phdi,wplm,
     *          j,m,nendyr,nbegyr,sx1,sx2,sx3,sx,indexj,indexm)
              GO TO 220
          ENDIF 
      ENDIF 
      PX2(j,m) = .897 * X2 + Z(j,m)/3.
      PX2(j,m) = AMIN1(PX2(j,m),0.)   
      IF (PX2(j,m).LE.-1.) THEN  
          IF (PX3(j,m).EQ.0.) THEN   
C IF NO EXISTING WET SPELL OR DROUGHT X2 BECOMES THE NEW X3 
              X(j,m)   = PX2(j,m) 
              PX3(j,m) = PX2(j,m) 
              PX2(j,m) = 0.  
              iass=2            
              CALL ASSIGN (iass,k8,ppr,px1,px2,px3,x,pdsi,phdi,wplm,
     *          j,m,nendyr,nbegyr,sx1,sx2,sx3,sx,indexj,indexm)
              GO TO 220
          ENDIF 
      ENDIF 
      IF (PX3(j,m).EQ.0.) THEN   
C NO ESTABLISHED DROUGHT (WET SPELL), BUT X3=0 SO EITHER X1 OR X2 USED AS X3   
          IF (PX1(j,m).EQ.0.) THEN   
              X(j,m) = PX2(j,m)   
              iass=2            
              CALL ASSIGN (iass,k8,ppr,px1,px2,px3,x,pdsi,phdi,wplm,
     *          j,m,nendyr,nbegyr,sx1,sx2,sx3,sx,indexj,indexm)
              GO TO 220
          ENDIF 
          IF (PX2(j,m).EQ.0) THEN
              X(j,m) = PX1(j,m)   
              iass=1   
              CALL ASSIGN (iass,k8,ppr,px1,px2,px3,x,pdsi,phdi,wplm,
     *          j,m,nendyr,nbegyr,sx1,sx2,sx3,sx,indexj,indexm)
              GO TO 220
          ENDIF 
      ENDIF 
C VALUES OF X1, X2, AND X3 ARE SAVED IN FILE 8. LATER X3 WILL BE THE 
C VALUE OF X (PDSI). THE ASSIGN SUBROUTINE BACKTRACKS THROUGH FILE 8 
C CHOOSING THE APPROPRIATE X1 OR X2 TO BE THAT MONTH'S X. 
      IF (K8.GT.40) STOP 'X STORE ARRAYS FULL'  
      SX1(K8) = PX1(j,m) 
      SX2(K8) = PX2(j,m) 
      SX3(K8) = PX3(j,m) 
      X(j,m)  = PX3(j,m) 
      K8  = K8 + 1
      k8max=k8  
      GO TO 220 
C PROB(END) RETURNS TO 0.  A POSSIBLE ABATEMENT HAS FIZZLED OUT,
C SO WE ACCEPT ALL STORED VALUES OF X3  
  210 PV  = 0.0 
      PX1(j,m) = 0.0 
      PX2(j,m) = 0.0 
      PPR(j,m) = 0.0 
      PX3(j,m) = .897 * X3 + Z(j,m)/3.
      X(j,m)   = PX3(j,m) 
      IF (K8.EQ.1) THEN 
        PDSI(j,m) = X(j,m)  
        PHDI(j,m) = PX3(j,m) 
        IF ( PX3(j,m).EQ. 0.0 ) PHDI(j,m)=X(j,m)
        CALL CASE (PPR(j,m),PX1(j,m),PX2(j,m),PX3(j,m),WPLM(j,m)) 
      ELSE  
        iass=3   
        CALL ASSIGN (iass,k8,ppr,px1,px2,px3,x,pdsi,phdi,wplm,
     *    j,m,nendyr,nbegyr,sx1,sx2,sx3,sx,indexj,indexm)
      ENDIF 
C SAVE THIS MONTH'S CALCULATED VARIABLES FOR USE WITH NEXT MONTHS DATA
  220 V   = PV 
      eff(j,m)=PV 
      PRO = PPR(j,m) 
      X1  = PX1(j,m) 
      X2  = PX2(j,m) 
      X3  = PX3(j,m) 
  230   continue
c Store these intermediate variables
2311   continue  
       if (iy.eq.nendyr) then
         call setmis(j,nendmo,eff,amsng)
         call setmis(j,nendmo,PPR,amsng)
         call setmis(j,nendmo,PX1,amsng)
         call setmis(j,nendmo,PX2,amsng)
         call setmis(j,nendmo,PX3,amsng)
       endif
       if(iy.ge.ioutbyr.and.iy.le.iouteyr)then
        write(iout(32),43)IST,IDV,iy,(eff(j,m),m=1,12)
        write(iout(33),43)IST,IDV,iy,(PPR(j,m),m=1,12)
        write(iout(34),43)IST,IDV,iy,(PX1(j,m),m=1,12)
        write(iout(35),43)IST,IDV,iy,(PX2(j,m),m=1,12)
        write(iout(36),43)IST,IDV,iy,(PX3(j,m),m=1,12)
       endif
C READ THE NEXT MONTH'S DATA   
  231 continue
C FINISH UP BY CLEARING OUT FILE 8 AND WRITING IT TO FILE 11
 240  do k8=1,k8max-1
      PDSI(indexj(k8),indexm(k8)) = X(indexj(k8),indexm(k8)) 
      PHDI(indexj(k8),indexm(k8)) = PX3(indexj(k8),indexm(k8)) 
      IF (PX3(indexj(k8),indexm(k8)).EQ. 0.0 )
     *   PHDI(indexj(K8),indexm(k8))=X(indexj(k8),indexm(k8))
      CALL CASE(PPR(indexj(k8),indexm(k8)),PX1(indexj(k8),indexm(k8)),
     *          PX2(indexj(k8),indexm(k8)),PX3(indexj(k8),indexm(k8)),
     *   WPLM(indexj(k8),indexm(k8))) 
      enddo 
c Call subroutine regid to get state or region character ID
      call regid(iSt,cSt)
C WRITE INDICES TO OUTPUT
      do j=1,nendyr-nbegyr+1
       iy=j+nbegyr-1
       write (75,118) cSt,IST,IDV,'05',iy,(PDSI(j,k),k=1,12),amsng
       write (76,118) cSt,IST,IDV,'06',iy,(PHDI(j,k),k=1,12),amsng
       write (77,118) cSt,IST,IDV,'07',iy,(Z(j,k),   k=1,12),amsng
       write (78,118) cSt,IST,IDV,'08',iy,(WPLM(j,k),k=1,12),amsng
  118  format (a2,I3.3,I2.2,a2,I4,13F7.2)
       if(iy.ge.ioutbyr.and.iy.le.iouteyr) then
        write(iout(11),43)IST,IDV,iy,(WPLM(j,m),m=1,12)
        write(iout(12),43)IST,IDV,iy,(PDSI(j,m),m=1,12)
        write(iout(13),43)IST,IDV,iy,(Z(j,m),m=1,12)
        write(iout(14),43)IST,IDV,iy,(PHDI(j,m),m=1,12)
        write(iout(15),43)IST,IDV,iy,(CP(j,m),m=1,12)
       endif
      enddo 
c Write current month's PHDI to current output file
      if ((cName(10)(13:17) .eq. 'div  ') .and. (iSt .le. 48)) then
       write(425,426) IST,IDV,PHDI(iouteyr-(nbegyr-1),iendmo)
 426   format(2I2.2,f7.2)
      endif
      goto 15
 16   continue
c Close all output files
      close (114)
      do 444,i=10,37
       if (i .le. 15) then
        close (iOut(i))
       else if (i.ge. 20) then
        close (iOut(i))
       endif
 444  continue
      do 76,i=75,78
       close (i)
 76   continue
      if (cName(10)(13:17) .eq. 'div  ') close (425)
      write(9,*) 'Subroutine pdinew executed'

      RETURN
      END   
c
c----------------- assign ---------------------------------------------
      SUBROUTINE ASSIGN (iass,k8,ppr,px1,px2,px3,x,pdsi,phdi,wplm,
     *          j,m,nendyr,nbegyr,sx1,sx2,sx3,sx,indexj,indexm)

      DIMENSION SX1(40), SX2(40), SX3(40), SX(96),
     *  pdsi(200,13),phdi(200,13),wplm(200,13),ppr(200,12),
     *  px1(200,12),px2(200,12),px3(200,12),x(200,12),
     *  indexj(40),indexm(40)
C   
C FIRST FINISH OFF FILE 8 WITH LATEST VALUES OF PX3, Z,X
C X=PX1 FOR I=1, PX2 FOR I=2, PX3,  FOR I=3 
      SX(K8) = X(j,m) 
      ISAVE  = iass
      IF (K8.EQ.1) THEN 
          PDSI(j,m) = X(j,m)  
          PHDI(j,m) = PX3(j,m) 
          IF ( PX3(j,m) .EQ. 0.0 ) PHDI(j,m)=X(j,m)
          CALL CASE (PPR(j,m),PX1(j,m),PX2(j,m),PX3(j,m),WPLM(j,m)) 
          RETURN
      ELSE  
      ENDIF 
      IF (Iass.EQ.3) THEN  
C USE ALL X3 VALUES
          DO 10 Mm = 1, K8-1 
          SX(Mm) = SX3(Mm)
   10     continue
          GO TO 70  
      ELSE  
C BACKTRACK THRU ARRAYS, STORING ASSIGNED X1 (OR X2) IN SX UNTIL IT IS ZERO,
C THEN SWITCHING TO THE OTHER UNTIL IT IS ZERO, ETC. 
          DO 60 Mm = K8-1, 1, -1 
          GO TO (20,40) ISAVE   
   20     IF (SX1(Mm).EQ.0) GO TO 50 
   30     ISAVE = 1 
          SX(Mm) = SX1(Mm)
          GO TO 60  
   40     IF (SX2(Mm).EQ.0) GO TO 30 
   50     ISAVE = 2 
          SX(Mm) = SX2(Mm)
   60     CONTINUE  
      ENDIF 
C PROPER ASSIGNMENTS TO ARRAY SX HAVE BEEN MADE OUTPUT THE MESS   
   70 DO 80 N = 1,K8
      PDSI(indexj(n),indexm(n)) = SX(N) 
      PHDI(indexj(n),indexm(n)) = PX3(indexj(n),indexm(n))
      IF(PX3(indexj(n),indexm(n)).EQ. 0.0)PHDI(indexj(n),indexm(n))
     *      =SX(N)
      CALL CASE (PPR(indexj(n),indexm(n)),PX1(indexj(n),indexm(n)),
     * PX2(indexj(n),indexm(n)),PX3(indexj(n),indexm(n)),
     * WPLM(indexj(n),indexm(n)))  
   80 CONTINUE  
      K8 = 1
      k8max=k8
      RETURN 
      END

c---------------- case -----------------------------------------------------
      SUBROUTINE CASE (PROB,X1,X2,X3,PALM)  
C   
C THIS SUBROUTINE SELECTS PRELIMINARY (OR NRT) PDSI FROM THE GIVEN X VALUES  
C AND THE PROBABILITY (PROB) OF ENDING EITHER A DROUGHT OR WET SPELL. 
C     X1   - INDEX FOR INCIPIENT WET SPELLS (ALWAYS POSITIVE)   
C     X2   - INDEX FOR INCIPIENT DRY SPELLS (ALWAYS NEGATIVE)   
C     X3   - SEVERITY INDEX FOR AN ESTABLISHED WET SPELL (POSITIVE) 
C            OR DROUGHT (NEGATIVE)  
C     PALM - THE SELECTED PDSI (EITHER PRELIMINARY OR FINAL)
C   This subroutine written and provided by CPC (Tom Heddinghaus
C   & Paul Sabol).
C   
      IF (X3.NE.0.) GO TO 10
C IF X3=0 THE INDEX IS NEAR NORMAL AND EITHER A DRY OR WET SPELL
C EXISTS.  CHOOSE THE LARGEST ABSOLUTE VALUE OF X1 OR X2.   
      PALM = X1 
      IF ( ABS(X2) .GT. ABS(X1) ) PALM=X2   
      GO TO 100 
C   
 10   CONTINUE  
      IF ( PROB.GT.0. .AND. PROB.LT.100. ) GO TO 20 
C A WEATHER SPELL IS ESTABLISHED AND PALM=X3 AND IS FINAL   
      PALM = X3 
      GO TO 100 
C   
 20   CONTINUE  
      PRO = PROB/100.   
      IF (X3 .GT. 0.) GO TO 30  
C TAKE THE WEIGHTED SUM OF X3 AND X1
      PALM = (1.-PRO)*X3 + PRO*X1   
      GO TO 100 
C   
 30   CONTINUE  
C TAKE THE WEIGHTED SUM OF X3 AND X2
      PALM = (1.-PRO)*X3 + PRO*X2   
 100  RETURN
      END   
c --------------------------------------------------------------------
       subroutine setmis(j,nendmo,d,amsng)
c   set to missing the values past the latest year-month with
c   data
       real d(200,12)
       if (nendmo.ge.12) return
       do 10 m=nendmo+1,12
       d(j,m)=amsng
10     continue
       return
       end
