      subroutine drotarea_update(LYR,LMON,IMULT)
c Subroutine to CALCULATE THE PERCENT AREA OF ANY REGION EXPERIENCING
c PALMER DROUGHT
C-WRITTEN 1-19-89 BY RICHARD HEIM JR. FOR THE NCDC HISTORICAL CLIMATE 
c Extensively modified, April 2016, Trevor Wallis
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,cYear*4,INFILE*30,
     :           cOfile(11)*9,cReg*2
       real AREA(48,10),rData(500,15),PMDIAR(200,500,15),dd(12,11),
     :      ARCAT(12,11),TPDSI(12),CATARA(50,200,12,11),dr(13)
       integer KOUNT(200),IST(500),IDIV(500),KREG(50)

      data cOfile /'paed41st.','pasd42st.','pamd43st.','pald44st.',
     :             'paid45st.','pann46st.','paiw47st.','pasw48st.',
     :             'pamw49st.','pavw50st.','paew51st.' /

      ANU = -9.99
 
       WRITE (9,*) 'Drought area processing started'
C READ DIVISIONAL AREA FILE
       open (3, file=cIndata//'divarea.', status='old')   
       DO 5 I=1,48
       READ (3,101) KK,(AREA(I,JJ),JJ=1,10)     
  5    continue     
101    FORMAT (I2,10F7.0) 
       close (3)
C READ PERIOD OF RECORD PMDI DIVISIONAL VALUES
       NYRS=0
       WRITE (9,*) 'reading p.o.r. divisional pmdi'
       WRITE (6,*) 'reading p.o.r. divisional pmdi'
       DO 10 IYR=1895,LYR
        call system('rm '//PATHDRDW//'pmdi.grep')
        write(cYear,'(i4.4)') IYR
        call system('grep '//cYear//' '//cOdata//
     :     'palmer/pmdi.index.div > '//PATHDRDW//'pmdi.grep')

        open(3,file=PATHDRDW//'pmdi.grep',status='old')
        do 100,i=1,500
         read(3,333,end=3333) (rData(i,iX),iX=1,14)
 333     format(f3.0,f2.0,5x,12f7.2)
         rData(i,15) = -9.99
 100    continue
 3333   close (3)
        KNT = i - 1
        IY=IYR-1894
        NYRS=IY
        KOUNT(IY)=KNT
        DO ID=1,KNT
        DO M=1,15
         PMDIAR(IY,ID,M)=rData(ID,M)
        enddo
        enddo
  10    continue
C READ CONTROL FILE
       open (2, file=PATHDRDW//'drotarea.update.control', status='old')
       NREGS=0
1      READ (2,200,END=299) IREG,INFILE
200    FORMAT (3X,I3,1X,A30)  
       NREGS=NREGS+1
       KREG(NREGS)=IREG
C READ FILE CONTAINING THE LIST OF STATE/DIVISIONS DEFINING THE REGION
       open(7,file=cIndata//INFILE,status='old')
       I=1  
       NDIVS=0    
       TAREA=0.0  
6      READ (7,201,END=99) IST(I),IDIV(I) 
201    FORMAT (I3,I2)   
       NDIVS=NDIVS + 1  
       TAREA=TAREA+AREA(IST(I),IDIV(I))   
       I=I+1
       GO TO 6 
 99    continue 
       CLOSE (7)
C COMPUTE THE P.O.R. % AREA VALUES FOR THIS REGION
       idd=0
       DO 20 IYR=1895,LYR 
        IY=IYR-1894    
C INITIALIZE ARRAYS
        DO II=1,12    
        DO JJ=1,11
         dd(ii,jj)=0.0    
         ARCAT(II,JJ)=0.0
        enddo
        enddo
        IENDMN=12     
        IF (IYR.EQ.LYR) IENDMN=LMON  
C Process into percent areas
       DO ID=1,KOUNT(IY)
       DO M=1,15
        rData(ID,M)=PMDIAR(IY,ID,M)
       enddo
       enddo
       DO 30 J=1,NDIVS  
       do 34 k=1,KOUNT(IY)
        is=INT(rData(k,1))
        id=INT(rData(k,2))
        if (is.eq.IST(J).and.id.eq.IDIV(J)) then
         DO 33 I=1,IENDMN
          TPDSI(I)=rData(k,I+2) 
          PP=TPDSI(I)
          IF (PP.GE.4.0) ARCAT(I,11)=ARCAT(I,11)+AREA(IST(J),IDIV(J))   
          IF (PP.GE.3.0.AND.PP.LT.4.0)
     *     ARCAT(I,10)= ARCAT(I,10)+AREA(IST(J),IDIV(J))  
          IF (PP.GE.2.0.AND.PP.LT.3.0)    
     *     ARCAT(I,9)= ARCAT(I,9)+AREA(IST(J),IDIV(J)) 
          IF (PP.GE.1.0.AND.PP.LT.2.0)    
     *     ARCAT(I,8)= ARCAT(I,8)+AREA(IST(J),IDIV(J)) 
          IF (PP.GE.0.5.AND.PP.LT.1.0)    
     *     ARCAT(I,7)= ARCAT(I,7)+AREA(IST(J),IDIV(J)) 
          IF (PP.GT.-0.5.AND.PP.LT.0.5)   
     *     ARCAT(I,6)= ARCAT(I,6)+AREA(IST(J),IDIV(J)) 
          IF (PP.GT.-1.0.AND.PP.LE.-0.5)  
     *     ARCAT(I,5)= ARCAT(I,5)+AREA(IST(J),IDIV(J)) 
          IF (PP.GT.-2.0.AND.PP.LE.-1.0)  
     *     ARCAT(I,4)= ARCAT(I,4)+AREA(IST(J),IDIV(J)) 
          IF (PP.GT.-3.0.AND.PP.LE.-2.0)  
     *     ARCAT(I,3)= ARCAT(I,3)+AREA(IST(J),IDIV(J)) 
          IF (PP.GT.-4.0.AND.PP.LE.-3.0)  
     *     ARCAT(I,2)=ARCAT(I,2)+AREA(IST(J),IDIV(J))     
          IF (PP.LE.-4.0) 
     *     ARCAT(I,1)=ARCAT(I,1)+AREA(IST(J),IDIV(J))
33       CONTINUE
         go to 30
        endif 
   34  continue

   30  continue  
C DETERMINE PERCENT OF TOTAL AREA FOR EACH CATEGORY
       DO M=1,IENDMN 
        DO N=1,11     
         dd(M,N)=100.*(ARCAT(M,N)/TAREA)
         CATARA(NREGS,IY,M,N)=dd(M,N)
        enddo
       enddo

 20   continue
C NEXT REGION
      GO TO 1
C DONE COMPUTING THE % AREA FOR P.O.R. FOR ALL REGIONS
299   CLOSE (2)
C PUT P.O.R. % AREAS FOR ALL REGIONS INTO FILES
       DO 800 ICAT=1,11
        open(95,file=cOdata//'palmer/'//cOfile(iCat),status='old')
        DO 801 IREG=1,NREGS
         call regid(kReg(iReg),cReg)
         DO 802 IYR=1895,LYR
          IY=IYR-1894
          write(95,95) cReg,kReg(iReg),'00',iCat+40,iYr,
     :      (CATARA(IREG,IY,M,ICAT),M=1,12),aNu
 95       format(a2,i3.3,a2,i2,i4,13f7.2)
 802     continue
 801    continue
        close (95)
 800   continue
C DONE
       write(9,*) 'Subroutine drotarea.update executed'
       RETURN
       END
c
