C weekmonthpalmerauto.f
c Program to compute and map weekly-monthly Palmer Index
C   WRITTEN 6-3-05 BY RICHARD HEIM JR. FOR NCDC/ScSD/CMB.
C Modified to operate on eddy with 1981-2010 mean values, August 2014,
c Trevor Wallis
C Major modification to clean up programs and operate in shared workspace,
c March 2016, Trevor Wallis
C23456789012345678901234567890123456789012345678901234567890123456789012

       COMMON /COMPATHS/ PATHCPC,PATHDRD,PATHDRDW,cMeanpath,PATHWDATA,
     :                   cOdata,cIndata
       CHARACTER PATHCPC*42,PATHDRD*34,PATHDRDW*40,cMeanpath*36,
     :           KOMMAND*200,INFILE*40,CDATE*24,FDATE*24,cPsdate*8,
     :           cCurrentDate*28,cIndata*38,MAIL*210,AMON(12)*3,
     :           cFtpsite*42,PATHWDATA*41,UPDATE*6,cOdata*38,cDay*3
       INTEGER NDAYS(12),JULIANL(1000) 
       LOGICAL FEXIST1,FEXIST2,DE_OUT

       DATA NDAYS/31,28,31,30,31,30,31,31,30,31,30,31/
       DATA AMON  / 'Jan','Feb','Mar','Apr','May','Jun',
     *              'Jul','Aug','Sep','Oct','Nov','Dec'  /
       PATHCPC=  '/workspace/drought/wkly_palmer/cpc1000day/'
       PATHWDATA='/workspace/drought/wkly_palmer/work_data/'
       PATHDRD=  '/workspace/nclimdiv/dev/work/data/'
       PATHDRDW= '/workspace/drought/wkly_palmer/programs/'
       cMeanpath='/workspace/1991-2020_means/clim_div/'
       cIndata = '/workspace/drought/wkly_palmer/indata/'
       cOdata =  '/workspace/drought/wkly_palmer/output/'
       cFtpsite = 'http://ftp.cpc.ncep.noaa.gov/htdocs/temp4/'

C23456789012345678901234567890123456789012345678901234567890123456789012
C Open processing log
       OPEN (9,FILE=PATHDRDW//'processing.log',STATUS='unknown')
c Make backup of pseudo.month.out
       call system('cp '//PATHWDATA//'pseudo.month.out '//
     :      PATHWDATA//'pseudo.month.out.bak')
C---GET DATE/TIME
       CALL SYSTEM ('rm -f '//PATHDRDW//'date.txt')
       CALL SYSTEM ('/bin/date > '//PATHDRDW//'date.txt')
       OPEN (29,FILE=PATHDRDW//'date.txt',STATUS='OLD')
       READ (29,2900) cCurrentDate
2900   FORMAT (A28)
       CLOSE (29)
C---OPEN OUTPUT MAIL FILE
       OPEN (29,FILE=PATHDRDW//'mailweekmon.',STATUS='UNKNOWN')
       WRITE (29,*) 'Output from weekmonthpalmerauto.exe started: ',
     :         cCurrentDate
C---CONVERT DATE-TIME
       cDay = cCurrentDate(1:3)
       DO I=1,12
         IF (cCurrentDate(5:7).EQ.AMON(I)) THEN
           JMON=I
           GO TO 10
         ENDIF
       ENDDO
2991   WRITE (29,*) '*** ERROR: Bad Date Format ',cCurrentDate
       go to 90
10     READ (cCurrentDate(9:10),'(I2)',ERR=2991) JDAY
       READ (cCurrentDate(25:28),'(I4)',ERR=2991) JYEAR
C---THE CPC FILES ARE TWO DAYS BEHIND.  DETERMINE DATE TWO DAYS AGO.
       IDAY=JDAY-2
       IMON=JMON
       IYEAR=JYEAR
c       IDAY=08
c       IMON=10
c       IYEAR=2022
       if (IDAY .lt. 1) then
        IMON = IMON - 1
        if (IMON .lt. 1) then
         IMON = 12
         IYEAR = IYEAR - 1
        endif
        if ((mod(IYEAR,4) .eq. 0) .and. (IMON .eq. 2)) then
         LDAY = 29
        else
         LDAY = NDAYS(IMON)
        endif
        IDAY = LDAY +IDAY
       endif
C---GET ust1000cd AND uspg1000cd FILES FROM CPC FTP SITE
c---delete old data files
ccc cf
ccc      call system('rm '//PATHCPC//'us*1000c*')
c---Get required files using WGET
ccc      call system('wget -q --directory-prefix='//PATHCPC//' '//
ccc     :       cFtpsite//'ust1000cd')
ccc      call system('wget -q --directory-prefix='//PATHCPC//' '//
ccc     :       cFtpsite//'uspg1000cd_new')
c     :       cFtpsite//'uspg1000cd')
ccc   cf
C---CREATE 1000 DAY FILE FROM nCLIMGRID DAILY
      call system(PATHDRDW//'Kdays.csh')
C---MAKE SURE FILES ust1000cd AND uspg1000cd ARE CURRENT
      INQUIRE (FILE=PATHCPC//'ust1000cd',EXIST=FEXIST1)
      INQUIRE (FILE=PATHCPC//'uspg1000cd_new',EXIST=FEXIST2)
***      INQUIRE (FILE=PATHCPC//'uspg1000cd',EXIST=FEXIST2)
      IF (FEXIST1.AND.FEXIST2) THEN
        WRITE (29,*) 'ust1000cd AND uspg1000cd_new EXIST'
***        WRITE (29,*) 'ust1000cd AND uspg1000cd EXIST'
      ELSE
        WRITE (29,*) 'ust1000cd AND uspg1000cd_new DO NOT EXIST'
***        WRITE (29,*) 'ust1000cd AND uspg1000cd DO NOT EXIST'
        go to 90
      ENDIF
c Calculate Julian day
      CALL JULIAN (NDAYS,IYEAR,IMON,IDAY,JULIANDAY)
      KYR=IYEAR-(IYEAR/100)*100
      KJULIAN=1000*KYR+JULIANDAY
      WRITE (29,2921) IYEAR,IMON,IDAY,JULIANDAY,KJULIAN
2921  FORMAT (5X,I4,2I3,I5,I10)
c Open and read 100 day precip data
       OPEN (1,FILE=PATHCPC//'uspg1000cd_new',STATUS='OLD')
***       OPEN (1,FILE=PATHCPC//'uspg1000cd',STATUS='OLD')
       II=0
       DO K=1,20
         READ (1,1001) (JULIANL(I+II),I=1,50)
1001     FORMAT (5X,50I5)
         II=II+50
       ENDDO
       CLOSE (1)
       DO K=1,1000
         IF (JULIANL(K).EQ.KJULIAN) GO TO 11
       ENDDO
       WRITE (29,*) 'ERROR: uspg1000cd_new NOT CURRENT'
***       WRITE (29,*) 'ERROR: uspg1000cd NOT CURRENT'
       write(9,*) 'ERROR: uspg1000cd_new NOT CURRENT'
***       write(9,*) 'ERROR: uspg1000cd NOT CURRENT'
       go to 90      
11     OPEN (1,FILE=PATHCPC//'ust1000cd',STATUS='OLD')
       II=0
       DO K=1,20
         READ (1,1001) (JULIANL(I+II),I=1,50)
         II=II+50
       ENDDO
       CLOSE (1)
       DO K=1,1000
         IF (JULIANL(K).EQ.KJULIAN) GO TO 12
       ENDDO
       WRITE (29,*) 'ERROR: ust1000cd NOT CURRENT'
       write(9,*) 'ERROR: ust1000cd NOT CURRENT'
       go to 90
12     CONTINUE
       WRITE (29,*) 'CPC files are current'
C---GENERATE PSEUDO-MONTHLY TEMPERATURE AND PRECIPITATION FOR THE 
C---DATE (WEEK) DESIRED AND CREATE pseudo.month.out FILE.
       CALL WEEKLYPSEUDOMONTHAUTO (NDAYS,IYEAR,IMON,IDAY)
       WRITE (29,*) 'Executed WEEKLYPSEUDOMONTHAUTO subroutine'
       WRITE (9,*) 'Executed WEEKLYPSEUDOMONTHAUTO subroutine'
C---COPY pseudo.month.out AS pseudo.month.out.YYYYMMDD AND COPY
C---IT TO working DIRECTORY
       write(cPsdate,'(i4.4,2i2.2)') IYEAR,IMON,IDAY
       call system('cp '//PATHWDATA//'pseudo.month.out '//
     :      PATHWDATA//'pseudo.month.out.'//cPsdate)
       WRITE (29,*) 'pseudo.month.out.'//cPsdate//' created'
       INFILE = 'pseudo.month.out.'//cPsdate
C---UPDATE pcpn01??. AND tmpc02??. FILES in working data directory
      if ((JMON .ne. 1) .and. (JDAY .ge. 6)) then
        call system('cp -f '//PATHDRD//'pcpn01??. '//PATHWDATA)
        call system('cp -f '//PATHDRD//'tmpc02??. '//PATHWDATA)
        WRITE (29,*) 'Current pcpn and tmpc files uploaded'
      ENDIF
C---CALL LOADWEEKLYDATA SUBROUTINE TO LOAD pseudo.month.out.YYYYMMDD 
       write (9,*) 'Calling LOADWEEKLYDATA subroutine'
       IMULT = 1
       UPDATE = ' '
       CALL LOADWEEKLYDATA (INFILE,IMULT,UPDATE)
       WRITE (29,*) 'LOADWEEKLYDATA subroutine executed'
       WRITE (9,*) 'LOADWEEKLYDATA subroutine executed'
       call system('rm '//PATHWDATA//INFILE)
C COMPUTE PALMER INDICES AND OTHER STATISTICS.
       WRITE (9,*) 'Executing Subroutine tkmean'
       call tkmean(IYEAR,IMON,IDAY,IMULT,UPDATE)
       WRITE (6,*) 'Subroutine tkmean executed'
c PALMER INDICES FOR DIVISIONS
       call system('cp '//cIndata//'soilconstdiv. '//
     :     PATHDRDW//'pdinew.soilconst')
       call system('cp '//cIndata//'pdinew.outputfiles.div '//
     :     PATHDRDW//'pdinew.outputfiles')
       write(9,*) 'Executing subroutine pdinew, divisional data'
       call pdinew(IYEAR,IMON,IDAY,IMULT,NDAYS)
       write(6,*) 'Subroutine pdinew executed, divisional data'
c PALMER INDICES FOR STATES AND REGIONS
       call system('cp '//cIndata//'soilconststate. '//
     :     PATHDRDW//'pdinew.soilconst')
       call system('cp '//cIndata//'pdinew.outputfiles.state '//
     :     PATHDRDW//'pdinew.outputfiles')
       write(9,*) 'Executing subroutine pdinew, states and regions'
       call pdinew(IYEAR,IMON,IDAY,IMULT,NDAYS)
       write(6,*) 'Subroutine pdinew executed, states and regions'
c PALMER INDICES FOR RIVER BASINS ETC
       call system('cp '//cIndata//'soilconstriver. '//
     :     PATHDRDW//'pdinew.soilconst')
       call system('cp '//cIndata//'pdinew.outputfiles.river '//
     :     PATHDRDW//'pdinew.outputfiles')
       write(9,*) 'Executing subroutine pdinew, river basins'
       call pdinew(IYEAR,IMON,IDAY,IMULT,NDAYS)
       write(6,*) 'Subroutine pdinew executed, river basins'
c REGIONAL PERCENT AREA AFFECTED BY DROUGHT OR WETNESS
       write(9,*) 'Executing subroutine drotarea_update'
       call drotarea_update(IYEAR,IMON,IMULT)
       write(6,*) 'Subroutine drotarea_update executed'
c Call divisional SPI subroutine
       write(9,*) 'Executing subroutine spinew for divisions'
       call spinew(IYEAR,1)
       write(6,*) 'Subroutine spinew executed for divisions'
       write(9,*) 'Executing subroutine spinew for states,etc.'
       call spinew(IYEAR,0)
       write(6,*) 'Subroutine spinew executed for states,etc.'
c Update prolog status and stored dates
       write(9,*) 'Calling update_prolog_status'
       call update_prolog_status()
       write(6,*) 'Subroutine update_prolog_status executed'
C RUN subroutine makemapsauto4 TO GENERATE MAPS
       write(9,*) 'Calling subroutine makemapsauto4 TO GENERATE MAPS'
       call makemapsauto4(IYEAR,IMON,IDAY,cDay)
       write(6,*) 'Subroutine makemapsauto4 has generated maps'
       write(9,*) 'Subroutine makemapsauto4 has generated maps'
       write(29,*) 'Subroutine makemapsauto4 has generated maps'
C RUN subroutine general_rank_auto to rank National temps
       write(9,*) 'Calling subroutine general_rank_auto'
       call general_rank_auto(IYEAR,IMON,IDAY)
       write(6,*) 'Subroutine general_rank_auto executed'
C SEND EMAIL
       write(29,*) 'Subroutine general_rank_auto executed'
 90    CLOSE (29)
       MAIL = 'mail -s "Weekly-Monthly Palmer output" '//
     :   ' Chris.Fenimore@noaa.gov < '//PATHDRDW//'mailweekmon.'
       CALL SYSTEM ('/bin/'//MAIL)
       write(6,*) 'Notification emailed'
       write(9,*) 'Notification emailed'

       STOP
       END

C23456789012345678901234567890123456789012345678901234567890123456789012
      SUBROUTINE getpor(nst,ndv,nel,dara)
C Gets temp or precip data
C
      COMMON /COMPATHS/ PATHCPC,PATHDRD,PATHDRDW,cMeanpath,PATHWDATA,
     :                  cOdata,cIndata
      CHARACTER PATHCPC*42,PATHDRD*34,PATHDRDW*40,cMeanpath*36,
     :          cElem(2)*4,cL*3,cSt*3,cDv*2,cEl*2,PATHWDATA*41,cCst*2,
     :          cOdata*38,cIndata*38
      real dara(200,14),amiss(2)

      DATA amiss / -9.99, -99.9 /
      DATA cElem /'pcpn','tmpc' /

c ---------------------------------------------------------------------
c Grep data
      if (ndv .eq. 0) then
       cL = 'st.'
      else
       cL = 'dv.'
      endif
c Set array tyo missing
      do 101,i=1,200
       dara(i,14) = -9999
       do 103,j=1,13
        dara(i,j) = amiss(nel)
 103   continue
 101  continue
c Get two character state or region abreviation, and form grep command
      call regid(nst,cCst)
      write(cSt,'(i3.3)') nst
      write(cDv,'(i2.2)') ndv
      write(cEl,'(i2.2)') nel
      call system('rm '//PATHDRDW//cElem(nel)//'.grep')
      call system('grep '//cCst//cSt//cDv//' '//PATHWDATA//
     :    cElem(nel)//cEl//cL//' > '//PATHDRDW//cElem(nel)//'.grep')
c ---------------------------------------------------------------------
c Open grepped file and read into array
      open(102,file=PATHDRDW//cElem(nel)//'.grep',status='old')
      do 100,i=1,200
       READ(102,16,end=18) dara(i,14),(dara(i,j),j=1,13)
 16    FORMAT(t10,f4.0,13f7.0)
 100  continue
 18   close (102)
      return
      end

C23456789012345678901234567890123456789012345678901234567890123456789012
       subroutine regid(nst,cCst)
c Subroutine to get two character state or region abbreviation
c ---------------------------------------------------------------------
c Type
      character cCst*2,cState(500)*2

      data cState /'AL','AZ','AR','CA','CO','CT','DE','FL','GA','ID',
     1             'IL','IN','IA','KS','KY','LA','ME','MD','MA','MI',
     2             'MN','MS','MO','MT','NE','NV','NH','NJ','NM','NY',
     3             'NC','ND','OH','OK','OR','PA','RI','SC','SD','TN',
     4             'TX','UT','VT','VA','WA','WV','WI','WY','**','AK',
     5             'HI','**','**','**','**','**','**','**','**','**',
     6             '**','**','**','**','**','PR','VI','**','**','**',
     7             '**','**','**','**','**','**','**','**','**','**',
     8             '**','**','**','**','**','**','**','**','**','**',
     9             'PI','**','**','**','**','**','**','**','**','**',
     *             'A1','A2','A3','A4','A5','A6','A7','A8','A9','AS',
     1             'AG','**','AN','AM','AU','**','**','**','**','AW',
     2             '**','**','**','**','**','**','**','**','**','**',
     3             '**','**','**','**','**','**','**','**','**','**',
     4             '**','**','**','**','**','**','**','**','**','**',
     5             'P1','P2','P3','P4','P5','P6','P7','P8','P9','PS',
     6             '**','**','**','**','**','**','**','**','**','**',
     7             '**','**','**','**','**','**','**','**','**','**',
     8             '**','**','**','**','**','**','**','**','**','**',
     9             '**','**','**','**','**','**','**','**','**','**',
     *             'Q1','Q2','Q3','Q4','Q5','Q6','Q7','Q8','Q9','QA',
     1             'QB','QC','QD','QE','QF','QG','QH','QI','**','QJ',
     2             '**','**','**','**','**','**','**','**','**','**',
     3             '**','**','**','**','**','**','**','**','**','**',
     4             '**','**','**','**','**','**','**','**','**','B1',
     5             '**','**','**','**','B2','B4','**','**','**','B3',
     6             'B5','B6','**','**','B7','**','**','**','**','**',
     7             '**','**','**','**','**','**','**','**','**','**',
     8             '**','**','**','**','**','**','**','**','**','**',
     9             '**','**','**','**','**','**','**','**','**','**',
     *             '**','**','**','**','**','**','**','**','**','**',
     1             '**','**','**','**','**','**','**','**','**','**',
     2             '**','**','**','**','**','**','**','**','**','**',
     3             '**','**','**','**','**','**','**','**','**','**',
     4             '**','**','**','**','**','**','**','**','**','E1',
     5             '**','**','**','**','**','E4','**','**','**','**',
     6             'E5','E6','**','**','E7','**','**','**','**','**',
     7             '**','**','**','**','**','**','**','**','**','**',
     8             '**','**','**','**','**','**','**','**','**','**',
     9             '**','**','**','**','**','**','**','**','**','**',
     *             '**','**','**','**','**','**','**','**','**','**',
     1             '**','**','**','**','**','**','**','**','**','**',
     2             '**','**','**','**','**','**','**','**','**','**',
     3             '**','**','**','**','**','**','**','**','**','**',
     4             '**','**','**','**','**','**','**','**','**','H1',
     5             '**','**','**','**','**','H4','**','**','**','**',
     6             'H5','H6','**','**','H7','**','**','**','**','**',
     7             '**','**','**','**','**','**','**','**','**','**',
     8             '**','**','**','**','**','**','**','**','**','**',
     9             '**','**','**','**','**','**','**','**','**','**'  /

      if ((nst .lt. 1) .or. (nst .gt. 500)) then
       write(6,*) nst,' is invalid regional code'
       stop
      endif
      cCst = cState(nst)
      if (cCst .eq. '**') then
       write(6,*) nst,' does not corespond to a real state or region'
       stop
      endif
 
      return
      end

C23456789012345678901234567890123456789012345678901234567890123456789012
       SUBROUTINE JULIAN (NDAYS,IYR,IMON,IDAY,JDAY)
C   Subroutine to get JULIAN DAY

       INTEGER NDAYS(12)

       JDAY=0
       if ( IMON .eq. 1) then
        JDAY = IDAY
       else
        DO 10 M=1,IMON-1
         if ((mod(IYR,4) .eq. 0) .and. (M .eq. 2)) then
          ND = 29
         else
          ND=NDAYS(M)
         endif
         JDAY = JDAY + ND
 10     continue
        JDAY = JDAY + IDAY
       endif 
       return
       END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
       SUBROUTINE WEEKLYPSEUDOMONTHAUTO (NDAYS,IYEAR,IMON,IDAY)
C   READ THE DAILY NORMAL DIVISIONAL FILES AND COMPUTE THE SUM OF THE
C   DAILY NORMAL PRECIPITATION AND TEMPERATURE FOR EACH DIVISION OVER
C   THE DESIRED TIME PERIOD (CONSECUTIVE DAYS WITHIN A SINGLE MONTH)
C
       COMMON /COMPATHS/ PATHCPC,PATHDRD,PATHDRDW,cMeanpath,PATHWDATA,
     :                   cOdata,cIndata
       COMMON /NRML1/ PCPNRML,TMPNRML 
       COMMON /CPC1/ PCP,TMP,ISTATE,IDIV
       CHARACTER PATHCPC*42,PATHDRD*34,PATHDRDW*40,cMeanpath*36,
     :           PATHWDATA*41,cOdata*38,cIndata*38
       REAL PCPNRML(48,10),TMPNRML(48,10),TMP(344),PCP(344)
       INTEGER ISTATE(344),IDIV(344),NDAYS(12)

       do 111,iSt=1,48
        do 112,iDv=1,10
         PCPNRML(iSt,iDv) = 0.
         TMPNRML(iSt,iDv) = 0.
 112    continue
 111   continue
       do 113,iDv=1,344
        TMP(iDv) = 0.
        PCP(iDv) = 0.
 113   continue
C---CONVERT THIS CALENDAR DATE TO A JULIAN DATE
       CALL JULIAN (NDAYS,IYEAR,IMON,IDAY,JULEND)
C---DETERMINE JULIAN DATE OF BEGINNING OF MONTH
       CALL JULIAN (NDAYS,IYEAR,IMON,1,JULBEG)
C---ADD DECADE AND UNIT PARTS OF YEAR TO JULIAN DATES
       JULBEG=1000*(IYEAR-(100*(IYEAR/100)))+JULBEG
       JULEND=1000*(IYEAR-(100*(IYEAR/100)))+JULEND
C---GET THE OBSERVED DATA SUMS FOR THE MONTH-TO-DATE
       CALL READDAT (JULBEG,JULEND)
       IDAYEND=NDAYS(IMON)
C---CHECK FOR LEAP YEAR
       IF ((IMON.EQ.2) .and. (MOD(IYEAR,4) .eq. 0)) IDAYEND=29
C---GET THE NORMAL DATA SUMS FOR THE REST OF THE MONTH
       IF (IDAY.LT.IDAYEND) THEN
         IDAYBEG=IDAY+1
         CALL READNRML (IMON,IDAYBEG,IDAYEND)
       ENDIF
C---ADD THE OBSERVED SUMS TO THE NORMAL SUMS & COMPUTE THE MEAN
C---TEMPERATURE FOR THE MONTH
       OPEN (UNIT=7,FILE=PATHWDATA//'pseudo.month.out',STATUS='UNKNOWN')
       WRITE (7,701) IYEAR,IMON,IDAY,JULEND,IDAY,IDAYEND-IDAY
701    FORMAT ('weekly ',I4,2I3.2,I7.5,3X,2I3)
       DO I=1,344
         PTOT=PCP(I)+PCPNRML(ISTATE(I),IDIV(I))
         TTOT=TMP(I)+TMPNRML(ISTATE(I),IDIV(I))
         TAVE=TTOT/FLOAT(IDAYEND)
         WRITE (7,702) ISTATE(I),IDIV(I),TAVE,PTOT,TMP(I),PCP(I),
     *     TMPNRML(ISTATE(I),IDIV(I)),PCPNRML(ISTATE(I),IDIV(I))
702      FORMAT (I3.3,I2.2,3(F7.1,F7.2))
       ENDDO
       CLOSE (7)
       RETURN
       END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
       SUBROUTINE READNRML (IMON,IBDAY,IEDAY)
C
C   READ THE DAILY NORMAL DIVISIONAL FILES AND COMPUTE THE SUM OF THE
C   DAILY NORMAL PRECIPITATION AND TEMPERATURE FOR EACH DIVISION OVER
C   THE DESIRED TIME PERIOD (CONSECUTIVE DAYS WITHIN A SINGLE MONTH)
C
       COMMON /NRML1/ PCPNRML,TMPNRML
       COMMON /COMPATHS/ PATHCPC,PATHDRD,PATHDRDW,cMeanpath,PATHWDATA,
     :                   cOdata,cIndata
       CHARACTER PATHCPC*42,PATHDRD*34,PATHDRDW*40,cMeanpath*36,
     :           PATHWDATA*41,cOdata*38,cIndata*38
       REAL PCPNRML(48,10),TMPNRML(48,10)
       INTEGER NDATA(31)
C---INITIALIZE
       DO IS=1,48
         DO ID=1,10
           PCPNRML(IS,ID)=0.0
           TMPNRML(IS,ID)=0.0
         ENDDO
       ENDDO 
C---READ PRECIPITATION DAILY NORMALS
       OPEN (UNIT=1,
     :  file=cMeanpath//'grdd.prcp.91-20.dly.means.dat',status='old')
1      READ (1,101,END=199) IST,IDV,MON,(NDATA(I),I=1,31)
101    format(2i2,2x,i2,31i5)
       IF (IST.GT.48) GO TO 1
       IF (MON.NE.IMON) GO TO 1
       IF (IBDAY.GT.IEDAY .OR. IEDAY.GT.31) THEN
         WRITE (6,*) 'DAY ERROR - READNRML PCP ',IBDAY,IEDAY
         STOP
       ENDIF
C---FOR FEB 29 (LEAP YEAR), USE NORMALS FROM FEB 28
       IF (IMON.EQ.2) NDATA(29)=NDATA(28)
C---COMPUTE SUM OF NORMAL PRECIPITATION OVER DESIRED PERIOD
       DO ID=IBDAY,IEDAY 
         PCPNRML(IST,IDV)=PCPNRML(IST,IDV)+FLOAT(NDATA(ID))/100.
       ENDDO
       GO TO 1
199    CLOSE (1)
C---READ TEMPERATURE DAILY NORMALS
       OPEN (UNIT=1,
     :  file=cMeanpath//'grdd.tave.91-20.dly.means.dat',status='old')
2      READ (1,101,END=299) IST,IDV,MON,(NDATA(I),I=1,31)
       IF (IST.GT.48) GO TO 2
       IF (MON.NE.IMON) GO TO 2
       IF (IBDAY.GT.IEDAY .OR. IEDAY.GT.31) THEN
         WRITE (6,*) 'DAY ERROR - READNRML TMP ',IBDAY,IEDAY
         STOP
       ENDIF
C---FOR FEB 29 (LEAP YEAR), USE NORMALS FROM FEB 28
       IF (IMON.EQ.2) NDATA(29)=NDATA(28)
C---COMPUTE SUM OF NORMAL TEMPERATURE OVER DESIRED PERIOD
       DO ID=IBDAY,IEDAY 
         TMPNRML(IST,IDV)=TMPNRML(IST,IDV)+FLOAT(NDATA(ID))/10.
       ENDDO
       GO TO 2
299    CLOSE (1)
       RETURN
       END 
C 
C23456789012345678901234567890123456789012345678901234567890123456789012
      SUBROUTINE READDAT (julbeg,julend)
c Subroutine to read CPC 1000 day files
C Original WRITTEN 2-13-2004 BY RICHARD HEIM JR
c Modified to operate on /monitoring-shared/workspace, Trevor Wallis
C January 2016
 
       COMMON /CPC1/ pcp,tmp,istate,idiv
       COMMON /COMPATHS/ PATHCPC,PATHDRD,PATHDRDW,cMeanpath,PATHWDATA,
     :                   cOdata,cIndata
       CHARACTER PATHCPC*42,PATHDRD*34,PATHDRDW*40,cMeanpath*36,
     :           PATHWDATA*41,cOdata*38,cIndata*38
       DIMENSION NDATA(1000,344),JUL(1000)
       integer istate(344),idiv(344)
       real tmp(344),pcp(344)

C---INITIALIZE/
       DO I=1,344
         tmp(I)=0.0
         pcp(I)=0.0
       ENDDO
c---read precipitation data
      open (unit=1,file=PATHCPC//'uspg1000cd_new',
***      open (unit=1,file=PATHCPC//'uspg1000cd',
     + blank='zero',status='unknown',access='sequential')
C   READ IN JULIAN DATES FROM 1000 DAY FILES
 906  FORMAT(5X,50I5)
      II=0
      DO 2 K=1,20
      READ(1,906) (JUL(I+II),I=1,50)
      II=II+50
   2  CONTINUE
C   LOOP THROUGH CLIMATIC DIVISIONS AND READ DATA
      DO 7 IDV=1,344            
      IS=1
      IE=50
      DO 4 I=1,20
      READ(1,102,END=16) IST,ICD,(NDATA(II,IDV),II=IS,IE)           
      IS=IS+50
      IE=IE+50
   4  CONTINUE
c---compute period total precipitation
       do 11 ii=1,1000
       if (JUL(ii).ge.julbeg .and. JUL(ii).le.julend) then
         pcp(IDV)=pcp(IDV)+(float(NDATA(ii,IDV))/100.)
       endif
11     continue
   7  CONTINUE
 102  FORMAT(1x,2I2,50I5)
16     close (1)  
c---read temperature data
      open (unit=1,file=PATHCPC//'ust1000cd',
     + blank='zero',status='unknown',access='sequential')
C   READ IN JULIAN DATES FROM 1000 DAY FILES                
      II=0
      DO 22 K=1,20
      READ(1,906) (JUL(I+II),I=1,50)
      II=II+50
  22  CONTINUE
C   LOOP THROUGH CLIMATIC DIVISIONS AND READ DATA
      DO 27 IDV=1,344
       ndays=0                                             
      IS=1
      IE=50
      DO 24 I=1,20
      READ(1,102,END=216) IST,ICD,(NDATA(II,IDV),II=IS,IE)
       istate(IDV)=IST
       idiv(IDV)=ICD           
      IS=IS+50
      IE=IE+50
  24  CONTINUE
c---compute period average temperature (COMPUTE SUM OF TEMP)
       do 211 ii=1,1000
       if (JUL(ii).ge.julbeg .and. JUL(ii).le.julend) then
         tmp(IDV)=tmp(IDV)+(float(NDATA(ii,IDV))/10.)
       endif
211    continue
  27  CONTINUE
216    close (1)
       RETURN
       END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
      SUBROUTINE LOADWEEKLYDATA (INFILE,IMULT,UPDATE)
C Subroutine to load weekly data
C WRITTEN 11-19-98 BY RICHARD HEIM JR. FOR NCDC CMB
c Modified to run on eddy using nCLIMDIV data and 1981-2010 means,
c August 2014, Trevor Wallis
c Modified to run on /monitoring-shared/workspace, Trevor Wallis Feb
c 2016
C Completely rewritten to comply with gFORTRAN, Trevor Wallis, March
C 2017

       COMMON /COMPATHS/ PATHCPC,PATHDRD,PATHDRDW,cMeanpath,PATHWDATA,
     :                   cOdata,cIndata
       CHARACTER PATHCPC*42,PATHDRD*34,PATHDRDW*40,cMeanpath*36,
     :           PATHWDATA*41,cOdata*38,cIndata*38,cDst(48)*2,
     :           cSt*2
       LOGICAL DE_OUT
       CHARACTER INFILE*40,UPDATE*6,CDATE*24,FDATE*24,cIndv(2)*9
       REAL rT(48,10),rP(48,10),rMiss(2),
     :      rDVar(48,10,1895:2050,13),rDNorm(48,10,12,13),rVar(13)

       data cIndv /'pcpn01dv.','tmpc02dv.'/
       data rMiss /-9.99,-99.9 /

       WRITE (9,*) 'LOADWEEKLYDATA processing started'
C---OPEN INPUT FILE AND READ HEADER RECORD
       OPEN (2,FILE=PATHWDATA//INFILE,status='old')
       READ (2,201) UPDATE,IYR,MON,IDAY
       write(6,201) UPDATE,IYR,MON,IDAY
201    FORMAT (A6,1X,I4,1X,I2,1X,I2)
       write (9,*) INFILE,' ',IYR,MON,IDAY
       IF (MON.LT.1 .OR. MON.GT.13) THEN
         WRITE (6,*) 'LOADWEEKLYDATA -- MONTH ERROR'
         STOP
       ENDIF
C---UPDATE payearmonweek. 
       OPEN (7,FILE=PATHDRDW//'payearmonweek.')
       READ (7,*) IYROLD,MONOLD,IDAYOLD
       CLOSE (7)
C---HAVE WE SKIPPED A MONTH OR YEAR?
       IF ((IYR.EQ.IYROLD .AND. MON.GT.MONOLD+1) .OR.
     *     (IYR.EQ.IYROLD+1 .AND. (MON.GT.1 .OR. MONOLD.LT.12)) .OR.
     *     (IYR.GT.IYROLD+1)) THEN
         WRITE (6,*) 'LOADWEEKLYDATA -- YEAR-MONTH ERROR'
         STOP
       ENDIF
C---UPDATE payearmonweek. FILE (USE TEMPORARY FILE FIRST)
       IF (UPDATE.EQ.'weekly' .OR. UPDATE.EQ.'WEEKLY') THEN
C---SET NEW-YEAR MULTIPLIER (1 = OLD YEAR, -1 = NEW YEAR)
        write (6,*) 'iyr, iyrold: ', iyr, iyrold
        write (6,*) 'mon, monold: ', mon, monold
        write (6,*) 'iday, idayold: ', iday, idayold
        IMULT=1
         IF (IYR.EQ.IYROLD+1) IMULT=-1
         IF ((IYR.EQ.IYROLD .AND. MON.EQ.MONOLD+1) .OR.
     *       (IYR.EQ.IYROLD .AND. MON.EQ.MONOLD .AND. 
     *             IDAY.GT.IDAYOLD) .OR.
     *       (IYR.EQ.IYROLD+1 .AND. MON.EQ.1 .AND. MONOLD.EQ.12)) THEN
           OPEN (7,FILE=PATHDRDW//'payearmon.temp',status='unknown')
           WRITE (7,701) IYR,MON,IDAY
701        FORMAT (I4,',',I2.2,',',I2.2)
           CLOSE (7)
         ELSE
           WRITE (9,953) IYR,MON,IYROLD,MONOLD,UPDATE
953        FORMAT ('LOADWEEKLYDATA ERROR: YEAR-MONTH ',2(I5,I3),1X,A12)
           WRITE (6,953) IYR,MON,IYROLD,MONOLD,UPDATE
ccc           STOP
         ENDIF
       ENDIF
c Read pseudo temps and precip from infile
c Initialize arrays
       do 100,iS=1,48
        do 101,iD=1,10
         rT(iS,iD) = -99.9
         rP(iS,iD) = -9.99
 101    continue
 100   continue
c Read data into arrays
       iCheck = 0
       do 110,i=1,1000
        read(2,202,END=299) iSt,iDv,rTemp,rPcp
202     FORMAT (I3,I2,F7.1,F7.2)
        rT(iSt,iDv) = rTemp
        rP(iSt,iDv) = rPcp
C---CHECK FOR DOB MISSING CODES
        IF (rT(iSt,iDv).lt.-90) rT(iSt,iDv) = -99.9
        IF (rP(iSt,iDv).lt.0) rP(iSt,iDv) = -9.99
C---CHECK THE MONTHLY TEMPERATURE AND PRECIPITATION FOR MISSING OR
C---INVALID VALUES
        IF ((rTemp.LE.-25) .or. (rTemp.ge.125)) then
         iCheck = iCheck + 1
         WRITE (9,903) IST,IDV,rTemp
903      FORMAT (10X,2I3,F7.1)
        endif
        IF ((rPcp.LT.0)  .or. (rPcp.GE.90)) then
         iCheck = iCheck + 1
         WRITE (9,904) IST,IDV,rPcp
904      FORMAT (10X,2I3,7X,F7.2)
        ENDIF
 110   continue
299    CLOSE (2)
       WRITE (9,*) 'FINISHED READING INPUT FILE'
C---ARE THERE ANY MISSING DATA FROM THE LOWER 48 STATES?
       IF (ICHECK.NE.0) THEN
         WRITE (9,910) IYR,MON,ICHECK
         WRITE (6,910) IYR,MON,ICHECK
910      FORMAT ('LOADWEEKLYDATA -- ',I4,I3,' HAS ',I4,
     *       ' ERRONEOUS VALUES. Check input data')
         STOP
       ENDIF
C---UPDATE update.status FILE
       OPEN (7,FILE=PATHDRDW//'update.status',STATUS='UNKNOWN')
       WRITE (7,702) IYR,MON,UPDATE,IDAY
702    FORMAT (I4,1X,I2,1X,A6,7X,I2)
       CLOSE (7)
       IF (UPDATE.EQ.'weekly' .OR. UPDATE.EQ.'WEEKLY') THEN
         ITYPE=1
       ELSE
         WRITE (9,*) 'UPDATE ERROR ',UPDATE
         STOP
       ENDIF
      WRITE (9,*) 'update.status updated',IYR,MON,UPDATE,IDAY
C---CALL SUBROUTINE PRLG1
       CALL PRLG1 (IYR,MON,ITYPE,4)
C---WRITE DIVISIONAL PRECIP AND TEMP TO DATA BASE
       DO 10 J=1,2
c Open POR infile and read into arrays
       open(107,file=PATHWDATA//cIndv(j),status='old')
c Initialize arrays
       do 120,iS=1,48
        cDst(iS) = '  '
        do 121,iD=1,10
         do 122,iM=1,13
          do 123,iY=1895,2050
           rDVar(iS,iD,iY,iM) = rMiss(j)
 123      continue
          do 124,iX=1,12
           rDNorm(iS,iD,iX,iM) = rMiss(j)
 124      continue
 122     continue
 121    continue
 120   continue
c Read input file
       do 130,iZ=1,100000
        if (j .eq. 1) then
         read(107,25,end=1000) cSt,iSt,iDv,iY,(rVar(iM),iM=1,13)
 25      format(a2,i3,i2,2x,i4,13f7.2)
        else
         read(107,26,end=1000) cSt,iSt,iDv,iY,(rVar(iM),iM=1,13)
 26      format(a2,i3,i2,2x,i4,13f7.1)
        endif
        if ((iSt .ge. 1) .and. (iSt .le. 48)) then
         cDst(iSt) = cSt
        else
         go to 130
        endif
c Sort Normal values
        if (iY .le. 9) then
         do 131,iM=1,13
          rDNorm(iSt,iDv,iY+1,iM) = rVar(iM)
 131     continue
        else if ((iY .eq. 31) .or. (iY .eq. 32)) then
         do 132,iM=1,13
          rDNorm(iSt,iDv,iY-20,iM) = rVar(iM)
 132     continue
        else if ((iY .ge. 1895) .and. (iY .le. 2050)) then
         do 133,iM=1,13
          rDVar(iSt,iDv,iY,iM) = rVar(iM)
 133     continue
        endif
 130   continue
 1000  close (107)
c Replace current month values with pseudo data
       do 140,iS=1,48
        do 141,iD=1,10
         if (j .eq. 1) then
          if (rP(iS,iD) .eq. rMiss(j)) go to 141
          rDVar(iS,iD,iYr,MON) = rP(iS,iD)
         else
          if (rT(iS,iD) .eq. rMiss(j)) go to 141
          rDVar(iS,iD,iYr,MON) = rT(iS,iD)
         endif
 141    continue
 140   continue
c Reopen divisional data files and write POR and pseudo data to it
       open(107,file=PATHWDATA//cIndv(j),status='unknown')
       do 150,iS=1,48
        do 151,iD=1,10
         if (j .eq. 1) then
          if (rP(iS,iD) .eq. rMiss(j)) go to 151
          do 152,iX=1,10
           write(107,27) cDst(iS),iS,iD,j,iX-1,
     :                   (rDNorm(iS,iD,iX,iM),iM=1,13)
 27        format(a2,i3.3,2i2.2,i4.4,13f7.2)
 152      continue
          do 153,iX=1,2
           write(107,27) cDst(iS),iS,iD,j,iX+30,
     :                   (rDNorm(iS,iD,iX+10,iM),iM=1,13)
 153      continue
          do 154,iY=1895,iYr
           write(107,27) cDst(iS),iS,iD,j,iY,
     :                   (rDVar(iS,iD,iY,iM),iM=1,13)
 154      continue
         else
          if (rT(iS,iD) .eq. rMiss(j)) go to 151
          do 155,iX=1,10
           write(107,28) cDst(iS),iS,iD,j,iX-1,
     :                   (rDNorm(iS,iD,iX,iM),iM=1,13)
 28        format(a2,i3.3,2i2.2,i4.4,13f7.1)
 155      continue
          do 156,iX=1,2
           write(107,28) cDst(iS),iS,iD,j,iX+30,
     :                   (rDNorm(iS,iD,iX+10,iM),iM=1,13)
 156      continue
          do 157,iY=1895,iYr
           write(107,28) cDst(iS),iS,iD,j,iY,
     :                   (rDVar(iS,iD,iY,iM),iM=1,13)
 157      continue
         endif
 151    continue
 150   continue
       close (107)
       WRITE (9,*) 'FINISHED UPDATING DATA FILE, ELEMENT ',J
10     CONTINUE
C---DONE
       IF (ICHECK.NE.0) THEN
         WRITE (9,907) IYR,MON,ICHECK
         WRITE (*,907) IYR,MON,ICHECK
907      FORMAT (' LOADWEEKLYDATA -- ',I4,I3,' HAS ',I4,
     *           ' POTENTIALLY ERRONEOUS VALUES')
       ENDIF
       IF (UPDATE.EQ.'weekly' .OR. UPDATE.EQ.'WEEKLY') THEN
         CALL SYSTEM ('cp -f '//PATHDRDW//'payearmon.temp '//
     *      PATHDRDW//'payearmon.')
         CALL SYSTEM ('mv -f '//PATHDRDW//'payearmon.temp '//
     *      PATHDRDW//'payearmonweek.')
         WRITE (9,*) 'payearmonweek. UPDATED: ',IYR,MON
       ELSE
         WRITE (6,*) 'payearmonweek. NOT UPDATED' 
         WRITE (9,*) 'payearmonweek. NOT UPDATED'
       ENDIF
       WRITE (9,*) 'LOADWEEKLYDATA -- SUCCESSFUL COMPLETION'
       WRITE (6,*) 'LOADWEEKLYDATA -- SUCCESSFUL COMPLETION'
       RETURN
       END
C ----------------------------------------------------------------------
       SUBROUTINE PRLG1 (IYR,MON,ITYPE,ISTATUS)
C
C   UPDATES ONE MONTH OF THE prolog.status FILE.
C   WRITTEN 11-19-98 BY RICHARD HEIM JR. FOR NCDC CMB
C
       COMMON /COMPATHS/ PATHCPC,PATHDRD,PATHDRDW,cMeanpath,PATHWDATA,
     :                   cOdata,cIndata
      CHARACTER PATHCPC*42,PATHDRD*34,PATHDRDW*40,cMeanpath*36,
     :           PATHWDATA*41,cOdata*38,cIndata*38
       LOGICAL CP_OUT
       INTEGER IARRAY(12,4)
C---DELETE TEMPORARY FILE IF IT ALREADY EXISTS
       INQUIRE (FILE=PATHDRDW//'prol.temp',EXIST=CP_OUT)
       IF (CP_OUT) CALL SYSTEM ('rm -f '//PATHDRDW//'prol.temp')
       OPEN (8,FILE=PATHDRDW//'prolog.status',STATUS='OLD')
       OPEN (10,FILE=PATHDRDW//'prol.temp',STATUS='NEW')
C---UPDATE CURRENT YEAR-MONTH IN prolog.status FILE
       DO 10 I=1,200
       READ (8,801,END=899) IY,((IARRAY(M,K),K=1,4),M=1,12)
801    FORMAT (1X,I4,12(1X,4I1))
       IF (IY.EQ.IYR) THEN
         IARRAY(MON,ITYPE)=ISTATUS
       ENDIF
       WRITE (10,1001) IY,((IARRAY(M,K),K=1,4),M=1,12)
1001   FORMAT ('Y',I4,12(1X,4I1))
10     CONTINUE
       WRITE (9,*) 'SUBROUTINE prlg1 -- TOO MANY YEARS ERROR'
       STOP
899    CLOSE (8)
       CLOSE (10)
       CALL SYSTEM 
     *   ('mv '//PATHDRDW//'prol.temp '//PATHDRDW//'prolog.status')
       WRITE (9,*) IYR,MON,ITYPE,ISTATUS,'prolog.status UPDATED'
       RETURN
       END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
       SUBROUTINE GETYEAR (nyear,ndv,nel,numrec,dara)
c Subroutine to return a year's data for one element
c Major modification of original subroutine, Trevor Wallis, March 2016
c ----------------------------------------------------------------------
c Type
       COMMON /COMPATHS/ PATHCPC,PATHDRD,PATHDRDW,cMeanpath,PATHWDATA,
     :                   cOdata,cIndata
      CHARACTER PATHCPC*42,PATHDRD*34,PATHDRDW*40,cMeanpath*36,
     :          PATHWDATA*41,cInfile(12)*9,cEl*2,cYear*4,
     :           grep_call*200,cOdata*38,cIndata*38
      real dara(500,15)

      data cInfile /'pcpn01dv.','tmpc02dv.','phdi06dv.',
     :              'zndx07dv.','pmdi08dv.','sp0171dv.',
     :              'sp0272dv.','sp0373dv.',
     :              'sp0674dv.','sp0975dv.','sp1276dv.','sp2477dv.' /

c ----------------------------------------------------------------------
      if (nel .eq. 1) then
       amiss = -9.99
      else
       amiss = -99.99
      endif
      write(cEl,'(i2.2)') nel
      write(cYear,'(i4.4)') nyear
      if ((nel .eq. 1) .or. (nel .eq. 2)) then
       grep_call = 'grep '//cEl//cYear//' '//PATHWDATA//
     :   cInfile(nel)//' > '//PATHDRDW//'temp'//cEl//'.file'
      else if ((nel .ge. 6) .and. (nel .le. 8)) then
       grep_call = 'grep '//cEl//cYear//' '//cOdata//'palmer/'//
     :   cInfile(nel-3)//' > '//PATHDRDW//'temp'//cEl//'.file'
      else if ((nel .ge. 71) .and. (nel .le. 77)) then
       grep_call = 'grep '//cEl//cYear//' '//cOdata//'spi/'//
     :   cInfile(nel-65)//' > '//PATHDRDW//'temp'//cEl//'.file'
      else
       write(6,*) 'GETYEAR has invalid Element input ',nel
       stop
      endif
c Initialize array
      do 13 i=1,500
       do 12 j=1,13
        dara(i,j)=amiss
   12   continue
        dara(i,14)=9999.
   13   continue
c ensure unique temporary file
      call system('rm '//PATHDRDW//'temp'//cEl//'.file')
      call system(grep_call)
      numrec=0
      open (102,file=PATHDRDW//'temp'//cEl//'.file')
      do 17 i=1,500
       read(102,16,end=18)(dara(i,j),j=1,15)
  16   format(2x,f3.0,f2.0,6x,13f7.2)
       numrec=numrec+1
  17  continue
  18  close(102)

      return
      end

C----------------------------------------------------------------------       
      SUBROUTINE SRMEAN(NDIVS,SCODE,SC,DIVS,STMEAN,iM)
      INTEGER SC
      REAL STMEAN,DIVS(53,10,12)
      REAL COEF(10,53) !AREA WEIGHTING COEFFICIENTS FOR EACH DIVISION   
      INTEGER NDIVS(53),SCODE(53)

      DATA ((COEF(I,K),I=1,10),K=1,53)/
     *.0893,.0974,.1482,.0775,.0970,.1495,.2852,.0559,.0000,.0000,
     *.1168,.3492,.0712,.0418,.0879,.1286,.2045,.0000,.0000,.0000,
     *.0901,.1043,.1389,.1266,.0987,.1319,.0952,.1258,.0885,.0000,
     *.1358,.1706,.0422,.0667,.2087,.0884,.2876,.0000,.0000,.0000,
     *.2704,.3695,.0833,.2033,.0735,.0000,.0000,.0000,.0000,.0000,
     *.2024,.6491,.1485,.0000,.0000,.0000,.0000,.0000,.0000,.0000,
     *.2209,.7791,.0000,.0000,.0000,.0000,.0000,.0000,.0000,.0000,
     *.2077,.2256,.1344,.2016,.1871,.0436,.0000,.0000,.0000,.0000,
     *.0632,.1036,.0574,.1192,.1581,.0989,.1167,.1404,.1425,.0000,
     *.0532,.0322,.0194,.3485,.0758,.1019,.0565,.0591,.0993,.1541,
     *.1284,.1169,.0985,.1095,.1013,.1336,.1331,.0960,.0827,.0000,
     *.1151,.1071,.0987,.1144,.1685,.0717,.1385,.1094,.0766,.0000,
     *.1217,.1083,.1172,.1309,.1307,.1077,.0887,.0991,.0957,.0000,
     *.0951,.1086,.0767,.0973,.1094,.1065,.1442,.1364,.1258,.0000,
     *.2438,.2271,.2120,.3171,.0000,.0000,.0000,.0000,.0000,.0000,
     *.0811,.1260,.0907,.0815,.1608,.1061,.1563,.0705,.1270,.0000,
     *.5394,.3103,.1503,.0000,.0000,.0000,.0000,.0000,.0000,.0000,
     *.1210,.1194,.1057,.0914,.0665,.3390,.0899,.0671,.0000,.0000,
     *.1911,.5107,.2982,.0000,.0000,.0000,.0000,.0000,.0000,.0000,
     *.1708,.1193,.0830,.1136,.0519,.0805,.0741,.0789,.1121,.1158,
     *.1562,.1575,.1227,.1142,.1212,.1074,.0729,.0768,.0711,.0000,
     *.0635,.1004,.0774,.0891,.1295,.1117,.1322,.0922,.1284,.0756,
     *.2110,.2008,.1473,.1652,.2165,.0592,.0000,.0000,.0000,.0000,
     *.1659,.0959,.1828,.1505,.1184,.1517,.1348,.0000,.0000,.0000,
     *.1842,.2438,.0990,.0000,.0954,.1087,.1227,.0589,.0873,.0000,
     *.2765,.3282,.3067,.0886,.0000,.0000,.0000,.0000,.0000,.0000,
     *.3078,.6922,.0000,.0000,.0000,.0000,.0000,.0000,.0000,.0000,
     *.3672,.5591,.0737,.0000,.0000,.0000,.0000,.0000,.0000,.0000,
     *.1119,.1522,.1237,.1182,.0794,.0791,.1800,.1555,.0000,.0000,
     *.0867,.1828,.1422,.0366,.1286,.0551,.0620,.0595,.1473,.0992,
     *.1464,.0680,.1182,.1061,.1087,.1862,.1339,.1325,.0000,.0000,
     *.1327,.1003,.1240,.1234,.1039,.0794,.1151,.1140,.1072,.0000,
     *.1188,.0803,.0918,.0986,.1169,.0721,.0801,.1101,.1003,.1310,
     *.1152,.1169,.1211,.0858,.1470,.1022,.0990,.1205,.0923,.0000,
     *.0863,.0973,.0918,.0302,.0620,.0918,.3058,.1313,.1032,.0000,
     *.0721,.0424,.0967,.0611,.1070,.0890,.0970,.0793,.1810,.1744,
     *1.000,.0000,.0000,.0000,.0000,.0000,.0000,.0000,.0000,.0000,
     *.0199,.1472,.1081,.2113,.1324,.1505,.2306,.0000,.0000,.0000,
     *.1639,.1020,.0991,.0367,.1333,.1434,.1130,.0761,.1325,.0000,
     *.2604,.1648,.3499,.2249,.0000,.0000,.0000,.0000,.0000,.0000,
     *.1474,.0956,.1459,.1242,.1353,.1320,.0803,.0457,.0821,.0115,
     *.2578,.0162,.0660,.2064,.1326,.0674,.2536,.0000,.0000,.0000,
     *.5051,.3306,.1643,.0000,.0000,.0000,.0000,.0000,.0000,.0000,
     *.1851,.1899,.1811,.1283,.1274,.1882,.0000,.0000,.0000,.0000,
     *.0739,.0266,.0530,.1412,.0910,.1388,.1028,.1739,.1222,.0766,
     *.0496,.1823,.2331,.3022,.1100,.1228,.0000,.0000,.0000,.0000,
     *.1727,.1776,.1102,.1262,.0966,.0835,.0985,.0866,.0481,.0000,
     *.0402,.0525,.2078,.1162,.1164,.0391,.0732,.1224,.0768,.1554,
     *.0564,.0296,.0193,.0372,.0786,.1069,.1483,.3091,.2147,.0000,
     *.1666,.1666,.1666,.1666,.1666,.1666,.0000,.0000,.0000,.0000,
     1.0830,.1205,.1537,.1764,.1675,.2989,.0000,.0000,.0000,.0000,
     *.3333,.3333,.3333,.0000,.0000,.0000,.0000,.0000,.0000,.0000,
     *.3333,.3333,.3333,.0000,.0000,.0000,.0000,.0000,.0000,.0000/

      DO 9 I=1,53
      IF(SC.EQ.SCODE(I))THEN
        ISC=I
        GOTO 8
      ENDIF
9     CONTINUE
       write (9,*) 'BAD STATE CODE in SRMEAN'
8     CONTINUE
      NUMD=NDIVS(ISC)
      STMEAN=0.0
      DO 10 I=1,NUMD
      STMEAN=STMEAN+DIVS(ISC,I,iM)*COEF(I,ISC)
10    CONTINUE
      RETURN
      END

C----------------------------------------------------------------------       
      SUBROUTINE update_prolog_status()
c Subroutine to update prolog.status and payearmon.

       COMMON /COMPATHS/ PATHCPC,PATHDRD,PATHDRDW,cMeanpath,PATHWDATA,
     :                   cOdata,cIndata
       CHARACTER PATHCPC*42,PATHDRD*34,PATHDRDW*40,cMeanpath*36,
     :           PATHWDATA*41,cOdata*38,cIndata*38,AY*1
       LOGICAL DE_OUT
       INTEGER IARRAY(12,4)

       INQUIRE(FILE=PATHDRDW//'prol.temp',EXIST=DE_OUT)
       IF (DE_OUT) CALL SYSTEM ('rm '//PATHDRDW//'prol.temp')
C---READ UPDATE STATUS INFORMATION
       OPEN (2,FILE=PATHDRDW//'payearmon.',STATUS='OLD')
       READ (2,*) IYR,MON
       CLOSE (2)
C---UPDATE prolog.status FILE
       OPEN (10,FILE=PATHDRDW//'prol.temp',STATUS='NEW')
       OPEN (7,FILE=PATHDRDW//'prolog.status',STATUS='OLD')
       DO 10 I=1,200
       READ (7,701,END=799) AY,IY,((IARRAY(M,K),K=1,4),M=1,12)
701    FORMAT (A1,I4,12(1X,4I1))
       IF (IY.EQ.IYR-1 .OR. IY.EQ.IYR) THEN
         IEMON=12
         IF (IY.EQ.IYR) IEMON=MON
         DO 11 M=1,IEMON
         DO 12 K=4,1,-1
         IF (IARRAY(M,K).EQ.1) THEN
           IARRAY(M,K)=2
           GO TO 11
         ENDIF
12       CONTINUE
11       CONTINUE
       ENDIF
       WRITE (10,701) AY,IY,((IARRAY(M,K),K=1,4),M=1,12)
10     CONTINUE

       WRITE (9,*) 'update.prolog.status TOO MANY YEARS ERROR'
       WRITE (6,*) 'update.prolog.status TOO MANY YEARS ERROR'
       stop

799    CLOSE (7)
       CLOSE (10)
       CALL SYSTEM ('mv '//PATHDRDW//'prol.temp '//
     :                  PATHDRDW//'prolog.status')
       WRITE (9,*) 'prolog.status UPDATED through ',IYR,MON
       RETURN
       END
************************************************************************
      subroutine tkmean(LYR,LMON,LDAY,IMULT,STATUS)

      COMMON /COMPATHS/ PATHCPC,PATHDRD,PATHDRDW,cMeanpath,PATHWDATA,
     :                   cOdata,cIndata
      CHARACTER PATHCPC*42,PATHDRD*34,PATHDRDW*40,cMeanpath*36,
     :           PATHWDATA*41,cOdata*38,STATUS*6,cIndata*38,
     :          cInfile(2)*9,cEl*2,cYear*4,grep_call*200,cOutfile(2)*9,
     :          cSt*2,cSst(48)*2,cRst(20)*2,cQst(20)*2,cCst(16)*2,
     :          cIndivs*20,cJunk*1
      INTEGER SCODE(53),NDIVS(53),MANY(9),LIST1(11),LIST2(4),LIST3(7),
     :        LIST4(6),LIST5(5),LIST6(6),LIST7(4),LIST8(3),LIST9(2)
      REAL REGION(12,9),STATES(12,53),SM(13),DIVS(53,10,12),TOTL(13),
     :     NATL(9),REG1(11),REG2(4),REG3(7),REG4(6),REG5(5),
     :     REG6(6),REG7(4),REG8(3),REG9(2),DVYR(500,15),rMiss(2),
     :     rSVar(48,1895:2050,13),rRVar(20,1895:2050,13),AREA(48,10),
     :     rQVar(20,1895:2050,13),rCVar(16,1895:2050,13),
     :     rSNorm(48,12,13),rRNorm(20,12,13),
     :     rQNorm(20,12,13),rCNorm(16,12,13),rVar(13)

      DATA REG1 / .02752,.01130,.18251,.05812,.04537,.05112,.04306,
     $          .27243,.24910,.00667,.05280/
      DATA REG2 /.22098,.22854,.33003,.22045/
      DATA REG3 /.18169,.11691,.13013,.22449,.13279,.13609,.07790/
      DATA REG4 /.17576,.19944,.20051,.17952,.10576,.13900/
      DATA REG5 /.31307,.16432,.15035,.16393,.20833/
      DATA REG6 /.09335,.14461,.08530,.08388,.12291,.46995/
      DATA REG7 /.26819,.24544,.28644,.19993/
      DATA REG8 /.33593,.38991,.27416/
      DATA REG9 /.58943,.41057/
      DATA NATL/.06021,.08428,.10271,.09715,.15551,.18822,.14053,
     $          .08230,.08908/
      DATA SCODE /1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,
     *  19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,
     *  37,38,39,40,41,42,43,44,45,46,47,48,50,51,66,67,91/
      DATA NDIVS /8,7,9,7,5,3,2,7,9,10,9,9,9,9,4,9,3,8,3,10,9,10,6,7,
     $           9,4,2,3,8,10,8,9,10,9,9,10,1,7,9,4,10,7,3,6,10,6,9,10,
     $           9,6,6,3,3/
      DATA MANY /11,4,7,6,5,6,4,3,2/
      DATA LIST1 /6,7,17,18,19,27,28,30,36,37,43/
      DATA LIST2 /13,20,21,47/
      DATA LIST3 /11,12,15,23,33,40,46/
      DATA LIST4 /01,08,09,31,38,44/
      DATA LIST5 /24,25,32,39,48/
      DATA LIST6 /03,14,16,22,34,41/
      DATA LIST7 /02,05,29,42/
      DATA LIST8 /10,35,45/
      DATA LIST9 /04,26/
      data cInfile /'pcpn01dv.','tmpc02dv.' /
      data cOutfile /'pcpn01st.','tmpc02st.' /
      data rMiss /-9.99,-99.9 /
c Eexecuting subroutine tkmean
      WRITE (9,*) 'Executing subroutine tkmean'
      IENDST=48
      IYRBM1=LYR-1
      open (2,file=cOdata//'tkmean.out',status='unknown')
C---FOR EACH ELEMENT (1,2):
      DO 1001 iEl=1,2
       write(cEl,'(i2.2)') iEl
       IFMULT=iEl*IMULT
c Read the old state data file into arrays. First initialize the arrays
       do 110,iM=1,13
        do 111,iY=1895,2050
         do 112,iSt=1,48
          rSVar(iSt,iY,iM) = rMiss(iEl)
 112     continue
         do 113,iR=1,20
          rRVar(iR,iY,iM) = rMiss(iEl)
 113     continue
         do 114,iQ=1,20
          rQVar(iQ,iY,iM) = rMiss(iEl)
 114     continue
         do 115,iC=1,16
          rCVar(iC,iY,iM) = rMiss(iEl)
 115     continue
 111    continue
        do 116,iV=1,12
         do 117,iSt=1,48
          rSNorm(iSt,iV,iM) = rMiss(iEl)
 117     continue
         do 118,iR=1,20
          rRNorm(iR,iV,iM) = rMiss(iEl)
 118     continue
         do 119,iQ=1,20
          rQNorm(iQ,iV,iM) = rMiss(iEl)
 119     continue
         do 120,iC=1,16
          rCNorm(iC,iV,iM) = rMiss(iEl)
 120     continue
 116    continue
 110   continue
       do 121,iSt=1,48
        cSst(iSt) = '  '
 121   continue
       do 122,iR=1,20
        cRst(iR) = '  '
 122   continue
       do 123,iQ=1,20
        cQst(iQ) = '  '
 123   continue
       do 124,iC=1,16
        cCst(iC) = '  '
 124   continue
c Open and read file
       open(107,file=PATHWDATA//cOutfile(iEl),status='old')
       do 130,iZ=1,100000
        if (iEl .eq. 1) then
         read(107,25,end=1500) cSt,iSt,iY,(rVar(iM),iM=1,13)
 25      format(a2,i3,4x,i4,13f7.2)
        else
         read(107,26,end=1500) cSt,iSt,iY,(rVar(iM),iM=1,13)
 26      format(a2,i3,4x,i4,13f7.1)
        endif
c sort into appropriate arrays
        if ((iSt .ge. 1) .and. (iSt .le. 48)) then
         cSst(iSt) = cSt
         do 131,iM=1,13
          if (iY .le. 9) then
           rSNorm(iSt,iY+1,iM) = rVar(iM)
          else if ((iY .eq. 31) .or. (iY .eq. 32)) then
           rSNorm(iSt,iY-20,iM) = rVar(iM)
          else if ((iY .ge. 1895) .and. (iY .le. 2050)) then
           rSVar(iSt,iY,iM) = rVar(iM)
          endif
 131     continue
        else if ((iSt .ge. 101) .and. (iSt .le. 120)) then
         cRst(iSt-100) = cSt
         do 132,iM=1,13
          if (iY .le. 9) then
           rRNorm(iSt-100,iY+1,iM) = rVar(iM)
          else if ((iY .eq. 31) .or. (iY .eq. 32)) then
           rRNorm(iSt-100,iY-20,iM) = rVar(iM)
          else if ((iY .ge. 1895) .and. (iY .le. 2050)) then
           rRVar(iSt-100,iY,iM) = rVar(iM)
          endif
 132     continue
        else if ((iSt .ge. 201) .and. (iSt .le. 220)) then
         cQst(iSt-200) = cSt
         do 133,iM=1,13
          if (iY .le. 9) then
           rQNorm(iSt-200,iY+1,iM) = rVar(iM)
          else if ((iY .eq. 31) .or. (iY .eq. 32)) then
           rQNorm(iSt-200,iY-20,iM) = rVar(iM)
          else if ((iY .ge. 1895) .and. (iY .le. 2050)) then
           rQVar(iSt-200,iY,iM) = rVar(iM)
          endif
 133     continue
        else if ((iSt .ge. 250) .and. (iSt .le. 265)) then
         cCst(iSt-249) = cSt
         do 134,iM=1,13
          if (iY .le. 9) then
           rCNorm(iSt-249,iY+1,iM) = rVar(iM)
          else if ((iY .eq. 31) .or. (iY .eq. 32)) then
           rCNorm(iSt-249,iY-20,iM) = rVar(iM)
          else if ((iY .ge. 1895) .and. (iY .le. 2050)) then
           rCVar(iSt-249,iY,iM) = rVar(iM)
          endif
 134     continue
        else
         go to 130
        endif
 130   continue
 1500  close (107)
C---INITIALIZE MISSING DATA MONTHS TO -9.99 OR -99.9
       IF (IEl.EQ.1) RH=-9.99
       IF (IEl.EQ.2) RH=-99.9
C---RECOMPUTE THE STATEWIDE, REGIONAL, NATIONAL VALUES FROM
C---BEGINNING OF RECORD TO PRESENT EACH RUN;
      write(6,*) 'Processing POR, 1895 - ',LYR
      DO 1000 IYR=1895,LYR
       write(cYear,'(i4.4)') iYr
C---GET ALL DIVISIONS FOR THIS YEAR
       call system('rm '//PATHDRDW//'gt'//cEl//'.file')
       grep_call = '/bin/grep '//cEl//cYear//' '//PATHWDATA//
     :   cInfile(iEl)//' > '//PATHDRDW//'gt'//cEl//'.file'
c Initialize array
       do 13 i=1,500
        DVYR(i,1) = 999.
        DVYR(i,2) = 99.
        do 12 j=3,15
         DVYR(i,j) = RH
 12     continue
 13    continue
c Make grep call
       call system(grep_call)
c Open and read grepped data file
       open(122,file=PATHDRDW//'gt'//cEl//'.file',status='old')
       KDVYR = 0
       do 17 i=1,500
        read(122,16,end=18)(DVYR(i,j),j=1,15)
  16    format(2x,f3.0,f2.0,6x,13f7.2)
        KDVYR = KDVYR + 1
  17   continue
  18   close(122)
       IF (IYR.LT.LYR) THEN
          IEMON=12
       ELSE IF (IYR.EQ.LYR) THEN
          IEMON=LMON
       ENDIF
       DO 837 IRHMS=1,12
       SM(IRHMS)=RH
       TOTL(IRHMS)=RH
837    continue
C---DO EACH STATE:  
      write(9,*) 'For each state. Year = ',iYr
      DO 1 ISC=1,IENDST
      IDV=NDIVS(ISC)
c For each division in state
      DO 10 ID=1,IDV
c Nebraska divn 4 does not exist
      IF(ISC.EQ.25.AND.ID.EQ.4)GOTO 10
       DO 14 IJK=1,KDVYR
       JSTA=INT(DVYR(IJK,1))
       JDIV=INT(DVYR(IJK,2))
       IF (JSTA.EQ.SCODE(ISC) .AND. JDIV.EQ.ID) THEN
         DO 20 IM=1,IEMON
         DIVS(ISC,ID,IM)=DVYR(IJK,IM+2)
20       continue
         GO TO 10
       ENDIF
14     CONTINUE
       WRITE (9,*) SCODE(ISC),ID,' NOT FOUND'
       WRITE (6,*) SCODE(ISC),ID,' NOT FOUND'
       STOP
10    CONTINUE
      DO 11 IM=1,IEMON
      CALL SRMEAN(NDIVS,SCODE,SCODE(ISC),DIVS,SM(IM),iM)
      STATES(IM,ISC)=SM(IM)
11    CONTINUE
      IF(IYR.GT.1980)WRITE(2,60)'SM:',iEl,SCODE(ISC),IDV,IYR,
     *     (SM(J),J=1,IEMON)
      if (iYr .eq. Lyr) rSVar(iSc,Lyr,Lmon) = States(Lmon,iSc)
60    FORMAT(1X,/, A10,I2,I3,I3,I5,/,1X,12F6.2)
1     CONTINUE
C---NOW DO REGIONS......
      DO 101 IREG=101,109
      IMANY=MANY(IREG-100)
      DO 401 IM=1,IEMON
      TOTL(IM)=0.0
401   continue
      DO 201 I=1,IMANY
      IF(IREG.EQ.101)THEN
        ISC=LIST1(I)
        XWT=REG1(I)
      ELSE IF(IREG.EQ.102)THEN
        ISC=LIST2(I)
        XWT=REG2(I)
      ELSE IF(IREG.EQ.103)THEN
        ISC=LIST3(I)
        XWT=REG3(I)
      ELSE IF(IREG.EQ.104)THEN
        ISC=LIST4(I)
        XWT=REG4(I)
      ELSE IF(IREG.EQ.105)THEN
        ISC=LIST5(I)
        XWT=REG5(I)
      ELSE IF(IREG.EQ.106)THEN
        ISC=LIST6(I)
        XWT=REG6(I)
      ELSE IF(IREG.EQ.107)THEN
        ISC=LIST7(I)
        XWT=REG7(I)
      ELSE IF(IREG.EQ.108)THEN
        ISC=LIST8(I)
        XWT=REG8(I)
      ELSE IF(IREG.EQ.109)THEN
        ISC=LIST9(I)
        XWT=REG9(I)
      ENDIF
      DO 501 IM=1,IEMON
      TOTL(IM)=TOTL(IM)+STATES(IM,ISC)*XWT
      REGION(IM,IREG-100)=TOTL(IM)
501   continue
201   CONTINUE
      IF(IYR.GT.1980)WRITE(2,61)'RE:',iEl,IREG,IYR,(TOTL(K),K=1,IEMON)
61     FORMAT(1X,/, A10,I2,I5,I6,/,1X,12F6.2)
      if (iYr .eq. LYR) then
       rRVar(IREG-100,LYR,LMON) = REGION(LMON,IREG-100)
      endif
101   CONTINUE
C---DO NATIONAL VALUE:
      DO 58 IM=1,IEMON
      TOTL(IM)=0.0
58    continue
      DO 59 IREG=101,109
      DO 160 IM=1,IEMON
       TOTL(IM)=TOTL(IM)+REGION(IM,IREG-100)*NATL(IREG-100)
160    continue
59    CONTINUE
      IF(IYR.GT.1980)WRITE(2,61)'NA:',iEl,110,IYR,(TOTL(K),K=1,IEMON)
      if (iYr .eq. LYR) then
       rRVar(10,LYR,LMON) = TOTL(LMON)
      endif
c ---------------------------------------------------------------------
c Now do additional regions. Only need this in current year
      if (iYr .ne. LYR) go to 1000
C READ DIVISIONAL AREA FILE
       open (3, file=cIndata//'divarea.', status='old')
       DO 5 I=1,48
       READ (3,199) KK,(AREA(I,JJ),JJ=1,10)
5      continue
199    FORMAT (I2,10F7.0)
       CLOSE (3)
c Open and step through wgtseas.control file
      open (4, file=PATHDRDW//'wgtseas.control', status='old')
      do 600,iW=1,32
*      do 600,iW=1,1
       read(4,599,end=499) cSt,iSt,cIndivs
 599   format(a2,i4,1x,a20)
c Open cIndivs
       open(3,file=cIndata//cIndivs,status='old')
c Read and reject header
       read(3,*) cJunk
c Read rest of file and process into an area weighted average
       n = 0
       rTotarea = 0
       rTotwtel = 0
       do 601,i=1,50
        read(3,37,end=1700) iSta,iDvn
 37     format(i3,i2)
        n = n + 1
        rTotarea = rTotarea + area(iSta,iDvn)
        rTotwtel = rTotwtel + divs(iSta,iDvn,LMON)*area(iSta,iDvn)
 601   continue
 1700  close (3)
       if (n .gt. 0) then
        rElmean = rTotwtel/rTotarea
       endif
c Write result to relavent regional array
       if ((iSt .gt. 100) .and. (iSt .le. 120)) then
        rRVar(iSt-100,LYR,LMON) = rElmean
       else if ((iSt .gt. 200) .and. (iSt .le. 220)) then
        rQVar(iSt-200,LYR,LMON) = rElmean
       else if ((iSt .ge. 250) .and. (iSt .le. 265)) then
        rCVar(iSt-249,LYR,LMON) = rElmean
       endif
 600  continue
 499  close (4)
c Go to next year
1000  continue
c rewrite state/region arrays to file
       open(107,file=PATHWDATA//cOutfile(iEl),status='unknown')
c states
       do 220,iSt =1,48
        if (cSst(iSt) .eq. '  ') go to 220
        do 221,iV=1,10
         if (iEl .eq. 1) then
          write(107,35) cSst(iSt),iSt,iEl,iV-1,
     :          (rSNorm(iSt,iV,iM),iM=1,13)
 35       format(a2,i3.3,2i4.4,13f7.2)
         else
          write(107,36) cSst(iSt),iSt,iEl,iV-1,
     :          (rSNorm(iSt,iV,iM),iM=1,13)
 36       format(a2,i3.3,2i4.4,13f7.1)
         endif
 221    continue
        do 222,iV=11,12
         if (iEl .eq. 1) then
          write(107,35) cSst(iSt),iSt,iEl,iV+20,
     :          (rSNorm(iSt,iV,iM),iM=1,13)
         else
          write(107,36) cSst(iSt),iSt,iEl,iV+20,
     :          (rSNorm(iSt,iV,iM),iM=1,13)
         endif
 222    continue
        do 223,iY=1895,LYR
         if (iEl .eq. 1) then
          write(107,35) cSst(iSt),iSt,iEl,iY,(rSVar(iSt,iY,iM),iM=1,13)
         else
          write(107,36) cSst(iSt),iSt,iEl,iY,(rSVar(iSt,iY,iM),iM=1,13)
         endif
 223    continue
 220   continue
c Regions
       do 320,iR=1,20
        if (cRst(iR) .eq. '  ') go to 320
        do 321,iV=1,10
         if (iEl .eq. 1) then
          write(107,35) cRst(iR),iR+100,iEl,iV-1,
     :          (rRNorm(iR,iV,iM),iM=1,13)
         else
          write(107,36) cRst(iR),iR+100,iEl,iV-1,
     :          (rRNorm(iR,iV,iM),iM=1,13)
         endif
 321    continue
        do 322,iV=11,12
         if (iEl .eq. 1) then
          write(107,35) cRst(iR),iR+100,iEl,iV+20,
     :          (rRNorm(iR,iV,iM),iM=1,13)
         else
          write(107,36) cRst(iR),iR+100,iEl,iV+20,
     :          (rRNorm(iR,iV,iM),iM=1,13)
         endif
 322    continue
        do 323,iY=1895,LYR
         if (iEl .eq. 1) then
          write(107,35) cRst(iR),iR+100,iEl,iY,
     :          (rRVar(iR,iY,iM),iM=1,13)
         else
          write(107,36) cRst(iR),iR+100,iEl,iY,
     :          (rRVar(iR,iY,iM),iM=1,13)
         endif
 323    continue
 320   continue
c River Basins ?
       do 420,iQ=1,20
        if (cQst(iQ) .eq. '  ') go to 420
        do 421,iV=1,10
         if (iEl .eq. 1) then
          write(107,35) cQst(iQ),iQ+200,iEl,iV-1,
     :          (rQNorm(iQ,iV,iM),iM=1,13)
         else
          write(107,36) cQst(iQ),iQ+200,iEl,iV-1,
     :          (rQNorm(iQ,iV,iM),iM=1,13)
         endif
 421    continue
        do 422,iV=11,12
         if (iEl .eq. 1) then
          write(107,35) cQst(iQ),iQ+200,iEl,iV+20,
     :          (rQNorm(iQ,iV,iM),iM=1,13)
         else
          write(107,36) cQst(iQ),iQ+200,iEl,iV+20,
     :          (rQNorm(iQ,iV,iM),iM=1,13)
         endif
 422    continue
        do 423,iY=1895,LYR
         if (iEl .eq. 1) then
          write(107,35) cQst(iQ),iQ+200,iEl,iY,
     :          (rQVar(iQ,iY,iM),iM=1,13)
         else
          write(107,36) cQst(iQ),iQ+200,iEl,iY,
     :          (rQVar(iQ,iY,iM),iM=1,13)
         endif
 423    continue
 420   continue
c Crop Regions
       do 520,iC=1,16
        if (cCst(iC) .eq. '  ') go to 520
        do 521,iV=1,10
         if (iEl .eq. 1) then
          write(107,35) cCst(iC),iC+249,iEl,iV-1,
     :          (rCNorm(iC,iV,iM),iM=1,13)
         else
          write(107,36) cCst(iC),iC+249,iEl,iV-1,
     :          (rCNorm(iC,iV,iM),iM=1,13)
         endif
 521    continue
        do 522,iV=11,12
         if (iEl .eq. 1) then
          write(107,35) cCst(iC),iC+249,iEl,iV+20,
     :          (rCNorm(iC,iV,iM),iM=1,13)
         else
          write(107,36) cCst(iC),iC+249,iEl,iV+20,
     :          (rCNorm(iC,iV,iM),iM=1,13)
         endif
 522    continue
        do 523,iY=1895,LYR
         if (iEl .eq. 1) then
          write(107,35) cCst(iC),iC+249,iEl,iY,
     :          (rCVar(iC,iY,iM),iM=1,13)
         else
          write(107,36) cCst(iC),iC+249,iEl,iY,
     :          (rCVar(iC,iY,iM),iM=1,13)
         endif
 523    continue
 520   continue
c ----------------------------------------------------------------------
c Close rewritten file
       close (107)

c Got to next element
1001  CONTINUE

      write (9,*) 'tkmean successfully updated. ',LYR,LMON
      write (6,*) 'tkmean successfully updated. ',LYR,LMON

      close (2)
      RETURN
      END

