C      ----------------------------------------------------------------
C      ----------------------------------------------------------------      
C                NOAA/Great Lakes Environmental Research Laboratory
C                            Ann Arbor, Michigan
C      -----------------------------------------------------------------
C      -----------------------------------------------------------------
C
C     Program Name      Midlakes
C     Language          FORTRAN
C     Compiler          Lahey 5.00
C
C     Source File       C:\coordmod\new\code\midlakes.for
C     Object File       C:\coordmod\new\code\midlakes.obj
C     Run Module        C:\coordmod\new\code\midlakes.exe
C
C     Subroutines       COEFF, MATRIX, STATS, STATSUB, OUTDATA,
C                       MONOUTLINE, STATOUTLINE, WRITEMONLINE,
C                       MONSTATSLINE, INITDATA, QREAD,
C                       NODAYS, CLRSCR, HALT
C      
C       This program is a hydrologic response model of the unregulated
C       portion of the Great Lakes (Michigan-Huron, St. Clair, and Erie,
C       and their connecting channels).  Quarter-monthly supplies are 
C       input either in the form of Net Basin Supplies or as hydrologic
C       components.  Input includes overlake precipitation, 
C       basin runoff, lake evaporation, diversions, groundwater, 
C       consumptive use, initial outflows and lake levels
C       for each lake and ice retardation rates and discharge equation
C       parameters for each of the connecting channels.  Outputs include
C       monthly mean lake levels, beginning-of-month levels, and
C       monthly outflows in the connecting channels.  Quarter-monthly 
C       Lake Ontario inflows are also calculated.  The Midlakes routing
C       model is comprised of continuity equations for each lake, solved
C       using a second-order finite-difference technique (Quinn, F.H.,1978.
C       Hydrologic Response Model of the North American Great Lakes,
C       Journal of Hydrology 37(3) pp. 295-307.).
C  
C       For further information contact:
C             Anne H. Clites
C             Great Lakes Environmental Research Laboratory
C             2205 Commonwealth Blvd.
C             Ann Arbor, Michigan  48105-1593
C             clites@glerl.noaa.gov
C    
C       Last revision: August 27, 1998
C

      PROGRAM MIDLAKES 
      
C  Variable declarations.   

      IMPLICIT NONE
      CHARACTER fnames(34)*40,title(3)*80,ext*4,type*4,run*3,dat*8,tim*8
      CHARACTER units*7
      INTEGER  qm, ierr, inc, i
      INTEGER  nyrtot, iyrn, iyr1, nyr, nyr1, nyear
      INTEGER  nfile, ndm, ndy, nruns
      INTEGER  mon, monb, mone, mon1, monn
      REAL mhmmbom(12,100),scmmbom(12,100),ermmbom(12,100)
      REAL mhmmflo(12,100),scmmflo(12,100),ermmflo(12,100)
      REAL mhmmlev(12,100),scmmlev(12,100),ermmlev(12,100)
      REAL mhlevsdv(12),mhbomsdv(12),mhflosdv(12),sclevsdv(12)
      REAL scbomsdv(12),scflosdv(12),erlevsdv(12),erbomsdv(12)
      REAL erflosdv(12),erbomsum(100)
      REAL mhlevsum(100),mhflosum(100),mhbomsum(100),sclevsum(100)
      REAL scflosum(100),scbomsum(100),erlevsum(100),erflosum(100)
      REAL amhlev(12),asclev(12),aerlev(12),amhbom(12),ascbom(12)
      REAL aerbom(12),amhflo(12),ascflo(12),aerflo(12)
      REAL zmhminle(12),zmhmaxle(12),zmhminbo(12),zmhmaxbo(12)
      REAL zmhminfl(12),zmhmaxfl(12),zscminle(12),zscmaxle(12)
      REAL zscminbo(12),zscmaxbo(12),zscminfl(12),zscmaxfl(12)
      REAL zerminle(12),zermaxle(12),zerminbo(12),zermaxbo(12)
      REAL zerminfl(12),zermaxfl(12),onflosum(100)
      REAL aonflo(12),onflosdv(12),zonmaxfl(12),zonminfl(12)
      REAL mhmnbs(12,100),scmnbs(12,100),ermnbs(12,100)
      REAL mhqnbs(4,12),scqnbs(4,12),erqnbs(4,12),niaqlev(4,12)
      REAL erqqflo(4,12,100)
      REAL onqinflo(4,12,100),onminflo(12,100)
      REAL miqpre(4,12), miqrun(4,12), miqevp(4,12), miqcon(4,12)
      REAL huqpre(4,12), huqrun(4,12), huqevp(4,12), huqcon(4,12)
      REAL scqpre(4,12), scqrun(4,12), scqevp(4,12), scqcon(4,12)
      REAL erqpre(4,12), erqrun(4,12), erqevp(4,12), erqcon(4,12)
      REAL miqgdw(4,12), huqgdw(4,12), scqgdw(4,12), erqgdw(4,12)
      REAL miqdiv(4,12), huqdiv(4,12), scqdiv(4,12), erqdiv(4,12)
      REAL mhqice(4,12), scqice(4,12), erqice(4,12), spqflo(4,12)
      REAL mhmflo,mhmnlev,mhlevmax,mhlevmin
      REAL mhmnflo,scmnflo,ermnflo
      REAL mhbommax,mhbommin,mhflomax,mhflomin
      REAL mhym,mhk,mhmlev,mhsar,misar,mhqmnlev,scqmnlev,erqmnlev
      REAL mha,mhb,mhwt,mhbom,mhqmnflo,scqmnflo,erqmnflo
      REAL mhlevann, sclevann, erlevann
      REAL mhbomann, scbomann, erbomann
      REAL mhfloann, scfloann, erfloann
      REAL mhz, scz, erz, niz, mhz2, scz2, erz2,deltaz4(4,12)
      REAL husar, scsar, ersar, erb, counta
      REAL scbom, erbom, scmnlev, ermnlev
      REAL dt, mhqflo, scqflo
      REAL sck, scym, sca, scb, scwt
      REAL erqflo, erk, erym, era, erwt
      REAL onqflo, scmlev, ermlev, scmflo, ermflo, onmflo
      REAL*8 c11, c12, c13, c21, c22, c23, c31, c32, c33
      REAL a1m1, a2m1, a3m1, b1m1, b2m1, b3m1
      REAL sdvmhlev, sdvmhbom, sdvmhflo, sdvsclev, sdvscbom, sdvscflo
      REAL sdverlev, sdverbom, sdverflo, sdvonflo
      REAL sclevmax, sclevmin, scbommax, scbommin, scflomax, scflomin
      REAL erlevmax, erlevmin, erbommax, erbommin, erflomax, erflomin
      REAL onfloann, onflomax, onflomin
      REAL*8 mhqwb,scqwb,erqwb
      REAL OUTFLOW
      REAL summhres,sumscres,sumerres,mhboq,scboq,erboq
      REAL summhft,sumscft,sumerft,mhft,scft,erft
      REAL*8 deltaz1,deltaz2,deltaz3


C  Common block statements.
      
      COMMON /BOP/      mhz,scz,erz,niz
      COMMON /COEFFS/   c11,c12,c13,c21,c22,c23,c31,c32,c33
      COMMON /EXP/      a1m1,a2m1,a3m1,b1m1,b2m1,b3m1
      COMMON /FILENAME/ fnames,nfile,nruns,run
      COMMON /HYDRO/    miqpre,huqpre,scqpre,erqpre,miqrun,huqrun,
     &                  scqrun,erqrun,miqdiv,huqdiv,scqdiv,erqdiv,
     &                  miqgdw,huqgdw,scqgdw,erqgdw,miqevp,huqevp,
     &                  scqevp,erqevp,miqcon,huqcon,scqcon,erqcon,
     &                  spqflo,mhqice,scqice,erqice,mhqnbs,scqnbs,
     &                  erqnbs,deltaz4,mhmnbs,scmnbs,ermnbs
      COMMON /LEVFLO/   mhmmlev,scmmlev,ermmlev,mhmmbom,scmmbom,ermmbom,
     &                  mhmmflo,scmmflo,ermmflo,onminflo,onqinflo,
     &                  mhbomsum,mhflosum,sclevsum,scbomsum,scflosum,
     &                  erlevsum,erbomsum,erflosum,mhlevsum,onflosum
      COMMON /MISC/     title,ext,type,dt,mon,nyr,inc,qm,units
      COMMON /NEWLEV/   mhz2,scz2,erz2
      COMMON /PERIODS/  iyr1,mon1,iyrn,monn
      COMMON /SA/       mhsar,misar,husar,scsar,ersar
      COMMON /STAGE/    mhk,sck,erk,mhym,scym,erym,mha,sca,era,
     &                  mhb,scb,erb,mhwt,scwt,erwt
      COMMON /STATS1/   amhbom,amhlev,amhflo,ascbom,asclev,ascflo,
     &                  aerbom,aerlev,aerflo,zmhmaxle,zmhminle,zmhmaxbo,
     &                  zmhminbo,zmhmaxfl,zmhminfl,zscmaxle,zscminle,
     &                  zscmaxbo,zscminbo,zscmaxfl,zscminfl,zermaxle,
     &                  zerminle,zermaxbo,zerminbo,zermaxfl,zerminfl
      COMMON /STATS2/   mhlevsdv,mhbomsdv,
     &                  mhflosdv,sclevsdv,scbomsdv,scflosdv,erlevsdv,
     &                  erbomsdv,erflosdv,sdvmhlev,sdvmhbom,sdvmhflo,
     &                  sdvsclev,sdvscbom,sdvscflo,sdverlev,sdverbom,
     &                  sdverflo,sdvonflo,mhlevann,mhbomann,mhfloann,
     &                  scbomann,scfloann,erlevann,erbomann,erfloann,
     &                  mhlevmax,mhlevmin,mhbommax,mhbommin,mhflomax,
     &                  mhflomin,sclevmax,sclevmin,scbommax,scbommin,
     &                  scflomax,scflomin,erlevmax,erlevmin,erbommax,
     &                  erbommin,erflomax,erflomin,aonflo,onfloann,
     &                  onflosdv,zonmaxfl,zonminfl,onflomax,onflomin,
     &                  sclevann
      COMMON /WATBAL/   mhqwb,scqwb,erqwb
      COMMON /DELTAS/ deltaz1,deltaz2,deltaz3



      WRITE (*,'(" Enter 3-letter file extension for this run:  ")')
      READ (5,2000) run

C  Date and Time stamp for 'details' file.

      CALL DATE(dat)
      CALL TIME(tim)

C  Open initial input data file and detailed quarter-monthly file.

      OPEN (UNIT=7, FILE='..\indata\initdata.'//run,STATUS='OLD',
     &      ACCESS='SEQUENTIAL', FORM='FORMATTED')
      OPEN (UNIT=35, FILE='..\outdata\details.'//run,STATUS='UNKNOWN',
     &      ACCESS='SEQUENTIAL', FORM='FORMATTED')

C  Read in initial data.

      WRITE (*,'(" Reading initdata file....")')
    1 CALL INITDATA
  
      WRITE(35,4000) 'details'//ext,dat,tim
      WRITE(35,4001)
      WRITE(35,4002) (title(i),i=1,3)
      WRITE(35,4003) inc
      WRITE(35,4004) mhk,mhym,mha,mhb,mhwt
      WRITE(35,4005) sck,scym,sca,scb,scwt
      WRITE(35,4006) erk,erym,era,erb,erwt

C  Calculate the number of scenario years for this run.
      
      nyrtot = iyrn - iyr1 + 1 
      nyr1=1

C  Initialize the annual statistics summation variables.

      mhbomann = 0.
      scbomann = 0.
      erbomann = 0.
      mhlevann = 0.
      sclevann = 0.
      erlevann = 0.
      mhfloann = 0.
      scfloann = 0.
      erfloann = 0.

C  Begin computations, one year at a time.
      
      DO nyr = nyr1, nyrtot
         nyear = iyr1 + nyr - nyr1

         PRINT*, nyear
         WRITE(35,4008)

C  Read in quarter-monthly net basin supplies for Mich-Hur, St. Clair
C  and Erie.  If 'type' = COMP, read in quarter-monthly components (precip,
C  runoff, and evaporation) and compute monthly net basin supplies.

         IF (type .NE. 'COMP') THEN

           CALL QREAD (nyear, fnames(1), mhqnbs)
           CALL QREAD (nyear, fnames(2), scqnbs)
           CALL QREAD (nyear, fnames(3), erqnbs)

           nfile = 3

         ELSE
           CALL QREAD (nyear, fnames(1), miqpre)
           CALL QREAD (nyear, fnames(2), huqpre)
           CALL QREAD (nyear, fnames(3), scqpre)
           CALL QREAD (nyear, fnames(4), erqpre)
           CALL QREAD (nyear, fnames(5), miqrun)
           CALL QREAD (nyear, fnames(6), huqrun)
           CALL QREAD (nyear, fnames(7), scqrun)
           CALL QREAD (nyear, fnames(8), erqrun)
           CALL QREAD (nyear, fnames(9), miqevp)
           CALL QREAD (nyear, fnames(10), huqevp)
           CALL QREAD (nyear, fnames(11), scqevp)
           CALL QREAD (nyear, fnames(12), erqevp)

           nfile = 12

C  Calculate quarter-monthly net basin supplies from components.
           
           DO mon = 1,12

              CALL NODAYS(nyear,mon,0,ndm,ndy,ierr)

              DO qm = 1,4

                 mhqnbs(qm,mon) = (miqpre(qm,mon)*misar +huqpre(qm,mon)*
     &              husar - miqevp(qm,mon)*misar -huqevp(qm,mon)*husar +
     &              miqrun(qm,mon)*misar + huqrun(qm,mon)*husar)
     &                  / (86400.*ndm/4.)

                 scqnbs(qm,mon) = (scqpre(qm,mon) - scqevp(qm,mon) + 
     &              scqrun(qm,mon))*scsar / (86400.*ndm/4.)

                 erqnbs(qm,mon) = (erqpre(qm,mon) - erqevp(qm,mon) + 
     &              erqrun(qm,mon))*ersar / (86400.*ndm/4.)

              ENDDO
           ENDDO

         ENDIF

C  Read in consumptive use and groundwater contribution values.
   
         CALL QREAD (nyear, fnames(nfile + 1), miqcon)
         CALL QREAD (nyear, fnames(nfile + 2), huqcon)
         CALL QREAD (nyear, fnames(nfile + 3), scqcon)
         CALL QREAD (nyear, fnames(nfile + 4), erqcon)

         CALL QREAD (nyear, fnames(nfile + 5), miqgdw)
         CALL QREAD (nyear, fnames(nfile + 6), huqgdw)
         CALL QREAD (nyear, fnames(nfile + 7), scqgdw)
         CALL QREAD (nyear, fnames(nfile + 8), erqgdw)

C  Read in diversions.

         CALL QREAD (nyear, fnames(nfile + 9), miqdiv)
         CALL QREAD (nyear, fnames(nfile + 10), huqdiv)
         CALL QREAD (nyear, fnames(nfile + 11), scqdiv)
         CALL QREAD (nyear, fnames(nfile + 12), erqdiv)

C  Read ice and weed retardation flows for connecting channels:
C  St. Clair R.(mhmice), Detroit R.(scmice), Niagara R.(ermice).
 
         CALL QREAD (nyear, fnames(nfile + 13), mhqice)
         CALL QREAD (nyear, fnames(nfile + 14), scqice)
         CALL QREAD (nyear, fnames(nfile + 15), erqice)

C  Read St. Marys River mean monthly flows. 

         CALL QREAD (nyear, fnames(nfile + 16), spqflo) 
  
C  If Niagara stage-discharge relationship is based on two gauges,
C  (erb greater than 0) read Chippawa-Grass Island Pool or Black Rock 
C  Lock beginning-of-month levels.

         IF (erb .GT. 0) THEN
            CALL QREAD (nyear, fnames(nfile + 17), niaqlev)
         ENDIF
 
         IF (ierr .EQ. 1) GO TO 900  

C  Initialize the yearly computational loop values.
      
         mhlevsum(nyr) = 0.  
         sclevsum(nyr) = 0.  
         erlevsum(nyr) = 0. 
         mhflosum(nyr) = 0. 
         scflosum(nyr) = 0.  
         erflosum(nyr) = 0.  
         onflosum(nyr) = 0.
         mhbomsum(nyr) = 0.
         scbomsum(nyr) = 0.
         erbomsum(nyr) = 0.

C  Initialize the counter for keeping track of how many months are
C  encountered in a single year, for statistics.
      
         counta = 0.

C  This is the monthly computational loop.

         monb = 1
         mone = 12
         IF (nyr .eq. nyr1)   monb=mon1
         IF (nyr .eq. nyrtot) mone=monn

         DO mon = monb, mone
 
            counta = counta + 1.

C  Store beginning-of-month lake levels: Michigan-Huron, St. Clair, Erie.

            mhmmbom(mon,nyr) = mhz 
            scmmbom(mon,nyr) = scz
            ermmbom(mon,nyr) = erz

C  Store beginning-of-month levels for annual statistics.

            mhbom = mhz
            scbom = scz
            erbom = erz

C  Initialize for mean monthly lake levels and flows.  Only half of 
C  the first and last timestep are used. These terms sum incremental
C  levels and flows for a month; a total of 4*inc increments.
      
            mhmnlev = .5 * mhz
            scmnlev = .5 * scz
            ermnlev = .5 * erz
            mhmnflo = .5 * OUTFLOW(mhz,scz,mhk,mhym,mha,mhb,mhwt)
            scmnflo = .5 * OUTFLOW(scz,erz,sck,scym,sca,scb,scwt)
            ermnflo = .5 * OUTFLOW(erz,niz,erk,erym,era,erb,erwt)

C  Call subroutine to determine number of days in the month.
C  Calculate the timestep in seconds.

            CALL NODAYS (nyear,mon,0,ndm,ndy,ierr)
            IF (ierr .EQ. 1) GO TO 900            
            dt = 86400.*ndm/4./inc

C  Quarter-monthly loop.

            DO qm = 1,4

               mhboq = mhz
               scboq = scz
               erboq = erz

C  Initialize quarter-month average levels and flows (half of first and
C  last segments).

               mhqmnlev = .5*mhboq
               scqmnlev = .5*scboq
               erqmnlev = .5*erboq
               mhqmnflo = .5*OUTFLOW(mhboq,scboq,mhk,mhym,mha,mhb,
     &                       mhwt)
               scqmnflo = .5*OUTFLOW(scboq,erboq,sck,scym,sca,scb,
     &                       scwt)
               erqmnflo = .5*OUTFLOW(erboq,niaqlev(qm,mon),erk,erym,
     &                       era,erb,erwt)
               
C  Compute quarter-monthly levels and outflows by increments.

               DO i = 1,inc
 
C  Subroutine Coeff uses levels of upstream and downstream lakes and 
C  stage-discharge parameters to calculate the 9 coefficients and water 
C  balance summation terms needed to solve the matrix.  Subroutine 
C  Matrix performs the matrix solution to determine change in levels.

                  CALL COEFF (i,mhz,scz,erz,niz,niaqlev)
                  CALL MATRIX (mhz,scz,erz)
                  
C  Update mean levels and beginning-of-increment levels.
C  Calculate quarter-monthly outflows by summing incremental flows.
            
		  mhqmnlev = mhqmnlev + mhz2
		  scqmnlev = scqmnlev + scz2
		  erqmnlev = erqmnlev + erz2
                  mhqflo = OUTFLOW(mhz2,scz2,mhk,mhym,mha,mhb,
     &              mhwt)
                  scqflo = OUTFLOW(scz2,erz2,sck,scym,sca,scb,
     &              scwt)
                  erqflo = OUTFLOW(erz2,niaqlev(qm,mon),erk,erym,
     &              era,erb,erwt)
                  mhqmnflo = mhqmnflo + mhqflo
                  scqmnflo = scqmnflo + scqflo
                  erqmnflo = erqmnflo + erqflo
                  mhmnlev = mhmnlev + mhz2
                  scmnlev = scmnlev + scz2
                  ermnlev = ermnlev + erz2
                  mhmnflo = mhmnflo + mhqflo
                  scmnflo = scmnflo + scqflo
                  ermnflo = ermnflo + erqflo
                  mhz = mhz2
                  scz = scz2
                  erz = erz2

               ENDDO

               niz=niaqlev(qm,mon)
  
C  End of incremental lake level computations loop.     
C  Compute quarter-monthly mean lake levels and flows (subtract half of last
C  segment for averaging). Subtract quarther-monthly ice/weed retardation from
C  quarter-monthly flows.

               mhqmnlev = (mhqmnlev - .5*mhz)/inc
               scqmnlev = (scqmnlev - .5*scz)/inc
               erqmnlev = (erqmnlev - .5*erz)/inc
               mhqmnflo = (mhqmnflo - .5*mhqflo)/inc - mhqice(qm,mon)
               scqmnflo = (scqmnflo - .5*scqflo)/inc - scqice(qm,mon)
               erqmnflo = (erqmnflo - .5*erqflo)/inc - erqice(qm,mon)

C  Note that erqmnflo, the Niagara outflow, does not include the Welland Canal
C  Save quarter-monthly flows for Erie needed for Ontario regulation model.

               erqqflo(qm,mon,nyr) = erqmnflo

C  Calculate quarter-monthly Ontario inflows for Ontario routing.  Ontario
C  inflows are Niagara outflows (erqqflo) plus the Welland Canal (erqdiv).
C  Diversion values are negative, so must be subtracted for Ontario inflow.      
               
               onqflo = erqqflo(qm,mon,nyr)-erqdiv(qm,mon)
               onqinflo(qm,mon,nyr) = onqflo

C  Print quarter-monthly results and inputs to 'details' file.

      WRITE(35,4007) nyear,mon,qm,mhqmnlev,mhqnbs(qm,mon),miqdiv(qm,mon)
     &  ,mhqice(qm,mon),mhqmnflo,scqmnlev,scqnbs(qm,mon),scqice(qm,mon),
     &  scqmnflo,erqmnlev,erqnbs(qm,mon),erqdiv(qm,mon),erqice(qm,mon),
     &  erqmnflo

            ENDDO

C  End of quarter-monthly loop.
C  Compute monthly mean lake levels and flows. Last segment divided as first was.
C  Average monthly ice/weed retardation subtracted from average monthly flow.
C  Since Ontario inflows are defined on a quarter-monthly basis (not incremental),
C  the monthly term is an average of quarter-monthlies.

            mhmlev = (mhmnlev - .5*mhz)/inc/4.
            scmlev = (scmnlev - .5*scz)/inc/4.                  
            ermlev = (ermnlev - .5*erz)/inc/4.
            mhmflo = (mhmnflo - .5*mhqflo)/inc/4. - (mhqice(1,mon)+
     &         mhqice(2,mon)+mhqice(3,mon)+mhqice(4,mon))/4.
            scmflo = (scmnflo - .5*scqflo)/inc/4. - (scqice(1,mon)+
     &         scqice(2,mon)+scqice(3,mon)+scqice(4,mon))/4.
            ermflo = (ermnflo - .5*erqflo)/inc/4. - (erqice(1,mon)+
     &         erqice(2,mon)+erqice(3,mon)+erqice(4,mon))/4.
            onmflo = (onqinflo(1,mon,nyr) + onqinflo(2,mon,nyr) +
     &         onqinflo(3,mon,nyr) + onqinflo(4,mon,nyr))/4.

C  Save monthly mean lake levels and flows.
      
            mhmmlev(mon,nyr) = mhmlev                        
            scmmlev(mon,nyr) = scmlev                       
            ermmlev(mon,nyr) = ermlev
            mhmmflo(mon,nyr) = mhmflo
            scmmflo(mon,nyr) = scmflo
            ermmflo(mon,nyr) = ermflo
            onminflo(mon,nyr) = onmflo

C  Save monthly average net basin supplies for summary file.

            mhmnbs(mon,nyr) = (mhqnbs(1,mon)+mhqnbs(2,mon)+mhqnbs(3,mon)
     &        +mhqnbs(4,mon))/4.

            scmnbs(mon,nyr) = (scqnbs(1,mon)+scqnbs(2,mon)+scqnbs(3,mon)
     &        +scqnbs(4,mon))/4.

            ermnbs(mon,nyr) = (erqnbs(1,mon)+erqnbs(2,mon)+erqnbs(3,mon)
     &        +erqnbs(4,mon))/4.

C  If units are English, convert from cfs to tcfs for printing to summary file.

            IF (units .NE. 'metric') THEN
               mhmnbs(mon,nyr) = mhmnbs(mon,nyr)/1000.
               scmnbs(mon,nyr) = scmnbs(mon,nyr)/1000.
               ermnbs(mon,nyr) = ermnbs(mon,nyr)/1000.
            ENDIF

C  Accumulate for annual mean levels and flows.
      
            mhlevsum(nyr) = mhlevsum(nyr) + mhmlev 
            sclevsum(nyr) = sclevsum(nyr) + scmlev 
            erlevsum(nyr) = erlevsum(nyr) + ermlev
            mhflosum(nyr) = mhflosum(nyr) + mhmflo 
            scflosum(nyr) = scflosum(nyr) + scmflo 
            erflosum(nyr) = erflosum(nyr) + ermflo
            onflosum(nyr) = onflosum(nyr) + onmflo
            mhbomsum(nyr) = mhbomsum(nyr) + mhbom
            scbomsum(nyr) = scbomsum(nyr) + scbom
            erbomsum(nyr) = erbomsum(nyr) + erbom

C  End of monthly loop.
           
            ENDDO

C  Advance to next year, determine annual mean lake levels and mean flows. 
  
         mhlevsum(nyr) = mhlevsum(nyr)/(counta*1.0)    
         mhbomsum(nyr) = mhbomsum(nyr)/(counta*1.0)    
         mhflosum(nyr) = mhflosum(nyr)/(counta*1.0)   
         sclevsum(nyr) = sclevsum(nyr)/(counta*1.0)     
         scbomsum(nyr) = scbomsum(nyr)/(counta*1.0)    
         scflosum(nyr) = scflosum(nyr)/(counta*1.0)     
         erlevsum(nyr) = erlevsum(nyr)/(counta*1.0)    
         erbomsum(nyr) = erbomsum(nyr)/(counta*1.0)     
         erflosum(nyr) = erflosum(nyr)/(counta*1.0)     
         onflosum(nyr) = onflosum(nyr)/(counta*1.0)     

C  End of yearly loop.

  100 ENDDO

C  Find out if another run is wanted using a new set of data.

      IF (nruns .GT. 1) GO TO 1

C  No new run wanted.  Call Subroutine Stats to perform statistical
C  calculations and then call Subroutine Outdata to open and print to
C  monthly output files for the three lakes. Each lake has 3 output files:
C  monthly mean levels (mhqlev.ext), monthly mean flows, (mhmflo.ext),
C  and beginning-of-quarter-month levels (mhblev.ext).

      CALL STATS
      CALL OUTDATA
         
C  Close all files

      CLOSE(7)
      CLOSE(8)
      CLOSE(9)
      CLOSE(10)
      CLOSE(11)
      CLOSE(12)
      CLOSE(13)
      CLOSE(21)
      CLOSE(22)
      CLOSE(23)
      CLOSE(31)
      CLOSE(32)
      CLOSE(33)
      CLOSE(34)
      CLOSE(35)

      STOP   

C  Error messages for data READ errors.

  900 WRITE(8,3000) nyear, mon, ndm
      STOP                     

C  Format statements.

 2000 FORMAT (A3)
 3000 FORMAT (5X,'NODAYS ERROR',3I5 )
 4000 FORMAT (A12,40X,A8,2x,A8)
 4001 FORMAT (/15X,'MIDLAKES - Quarter-Monthly Details file'/)
 4002 FORMAT (A80)
 4003 FORMAT (/'Time increments per month:',I3/)
 4004 FORMAT ('Stage-discharge equation parameters:'/20X,'-K-',7x,'-YM-'
     &,5x,'-A-',4X,'-B-',4X,'-WT-',/,'Michigan-Huron: ',F9.4,4X,F6.2,4X,
     &F3.1,4X,F3.1,4X,F3.1)
 4005 FORMAT ('St. Clair: ',5X,F9.4,4X,F6.2,4X,F3.1,4X,F3.1,4X,F3.1)
 4006 FORMAT ('Erie: ',10X,F9.4,4X,F6.2,4X,F3.1,4X,F3.1,4X,F3.1,4X)
 4007 FORMAT (I4,1X,I2,1X,I1,1X,F7.3,1X,F8.0,1X,F6.0,1X,F6.0,1X,F8.0,2X,
     &F7.3,1X,F8.0,1X,F6.0,1X,F8.0,2X,F7.3,1X,F8.0,1X,F6.0,1X,
     &F6.0,1X,F8.0)
 4008 FORMAT (/10X,'mhqlev',4x,'mhqnbs',1X,'mhqdiv',1X,'mhqice',2X,
     &'mhqflo',3X,'scqlev',4x,'scqnbs',1x,'scqice',2X,'scqflo',3X,
     &'erqlev',4X,'erqnbs',1X,'erqdiv',1X,'erqice',2X,'erqflo'/)

      END

C----------------------------------------------------------------------
      Function OUTFLOW(zu,zd,k,ym,a,b,wt)
C     -----------------------------------
C
C  Given constants k, ym, a, and b, and variables zu and zd, 
C  calculate channel outflow.  Terms include:
C         zu = upstream lake elevation
C         zd = downstream lake elevation
C         k = outflow equation coefficient
C         ym = outlet bottom elevation
C         a = depth exponent
C         b = fall exponent
C         wt = relative weight between upstream and downstream lake levels
C  The constants k, ym, a, and b are dictated by the stage-discharge 
C  relationship for each connecting channel.  This function was written by
C  Anne H. Clites, NOAA/GLERL, September, 1995.

      REAL zu,zd
      REAL k,ym,a,b,wt,wtdepth,fall,outflow

C     Error checking.

      IF (zd .GT. zu) THEN
         WRITE(*,'(1X,/)')
         WRITE(*,'(1X,"zd=",F6.2,"zu=",F6.2)') zd, zu
         CALL HALT("***zd is greater than zu")
      END IF
      IF (ym .GT. zu) THEN
         WRITE(*,'(/)')
         WRITE(*,'(1X,"ym=",F6.2,"zu=",F6.2)') ym, zu
         CALL HALT("***ym is greater than zu")
      END IF

      wtdepth = wt*zu + (1.-wt)*zd - ym
      fall = zu-zd
      outflow = k*(wtdepth)**a*(fall)**b
      RETURN
      END
C----------------------------------------------------------------------
      Subroutine COEFF(i,mhz,scz,erz,niz,niaqlev)
C     -------------------------------------
C
C     This subroutine calculates the 9 matrix coefficients for a specified
C     combination of gauges. When b=0, there is a single gauge; when b is 
C     greater than 0, there is more than one gauge involved.  This subroutine  
C     was written by Anne H. Clites, NOAA/GLERL, October, 1995.
C
      CHARACTER*80 title(3), fnames(34)*40, ext*4, run*3, units*7
      INTEGER qm
      REAL mhk,mhym,mhsar,misar
      REAL deltaz4(4,12)
      REAL niaqlev(4,12),mhqnbs(4,12),scqnbs(4,12),erqnbs(4,12)
      REAL miqpre(4,12),miqrun(4,12),miqevp(4,12),miqcon(4,12)
      REAL huqpre(4,12), huqrun(4,12), huqevp(4,12), huqcon(4,12)
      REAL scqpre(4,12), scqrun(4,12), scqevp(4,12), scqcon(4,12)
      REAL erqpre(4,12), erqrun(4,12), erqevp(4,12), erqcon(4,12)
      REAL miqgdw(4,12), huqgdw(4,12), scqgdw(4,12), erqgdw(4,12)
      REAL miqdiv(4,12), huqdiv(4,12), scqdiv(4,12), erqdiv(4,12)
      REAL mhqice(4,12), scqice(4,12), erqice(4,12), spqflo(4,12)
      REAL mhmnbs(12,100),scmnbs(12,100),ermnbs(12,100)
      REAL OUTFLOW,out1,out2,out3,out4,out5,out6
      REAL mha,sca,era,mhb,scb,erb,mhwt,scwt,erwt
      REAL a1m1,a2m1,a3m1,b1m1,b2m1,b3m1
      REAL*8 c11,c12,c13,c21,c22,c23,c31,c32,c33
      REAL*8 mhqwb,scqwb,erqwb
      REAL mhz,scz,erz,niz

C     Common block statements.

      COMMON /COEFFS/  c11,c12,c13,c21,c22,c23,c31,c32,c33
      COMMON /EXP/     a1m1,a2m1,a3m1,b1m1,b2m1,b3m1
      COMMON /FILENAME/ fnames,nfile,nruns,run
      COMMON /HYDRO/   miqpre,huqpre,scqpre,erqpre,miqrun,huqrun,
     &                 scqrun,erqrun,miqdiv,huqdiv,scqdiv,erqdiv,
     &                 miqgdw,huqgdw,scqgdw,erqgdw,miqevp,huqevp,
     &                 scqevp,erqevp,miqcon,huqcon,scqcon,erqcon,
     &                 spqflo,mhqice,scqice,erqice,mhqnbs,scqnbs,
     &                 erqnbs,deltaz4,mhmnbs,scmnbs,ermnbs
      COMMON /MISC/    title,ext,type,dt,mon,nyr,inc,qm,units
      COMMON /PERIODS/ iyr1,mon1,iyrn,monn            
      COMMON /SA/      mhsar,misar,husar,scsar,ersar
      COMMON /STAGE/   mhk,sck,erk,mhym,scym,erym,mha,sca,era,
     &                 mhb,scb,erb,mhwt,scwt,erwt
      COMMON /WATBAL/  mhqwb,scqwb,erqwb     

      a1m1=mha-1.
      a2m1=sca-1.
      a3m1=era-1.
      b1m1=mhb-1.
      b2m1=scb-1.
      b3m1=erb-1.
 
C  Calls to Function OUTFLOW needed in this subroutine.

      outmh=OUTFLOW(mhz,scz,mhk,mhym,mha,mhb,mhwt)  
      out1=OUTFLOW(mhz,scz,mhk,mhym,a1m1,mhb,mhwt)
      out2=OUTFLOW(mhz,scz,mhk,mhym,mha,b1m1,mhwt)
      outsc=OUTFLOW(scz,erz,sck,scym,sca,scb,scwt)
      out3=OUTFLOW(scz,erz,sck,scym,a2m1,scb,scwt)
      out4=OUTFLOW(scz,erz,sck,scym,sca,b2m1,scwt)
      IF(i.eq.1)THEN
        outer=OUTFLOW(erz,niz,erk,erym,era,erb,erwt) 
        out5=OUTFLOW(erz,niz,erk,erym,a3m1,erb,erwt)
        out6=OUTFLOW(erz,niz,erk,erym,era,b3m1,erwt)
      ELSE
        outer=OUTFLOW(erz,niaqlev(qm,mon),erk,erym,era,erb,erwt) 
        out5=OUTFLOW(erz,niaqlev(qm,mon),erk,erym,a3m1,erb,erwt)
        out6=OUTFLOW(erz,niaqlev(qm,mon),erk,erym,era,b3m1,erwt)
      END IF

C  Calculate matrix coefficients for this time step.

      c11 = .5*mha*mhwt*out1 + .5*mhb*out2 + mhsar/dt
      c12 = .5*mha*(1.-mhwt)*out1 - .5*mhb*out2  
      c13 = 0.
      c21 = -(.5*mha*mhwt*out1+.5*mhb*out2)
      c22 = .5*sca*scwt*out3 + .5*scb*out4 - .5*mha*(1.-mhwt)*out1 -
     &      .5*mhb*out2 + scsar/dt 
      c23 = -(-.5*sca*(1.-scwt)*out3 + .5*scb*out4)
      c31 = 0.
      c32 = -(.5*sca*scwt*out3+.5*scb*out4) 
      c33 = .5*era*erwt*out5 + .5*erb*out6 - 
     &      .5*sca*(1.-scwt)*out3 + .5*scb*out4 + ersar/dt 

C  Calculate deltaz4, the change in levels of the Chippewa-
C  Grass Island Pool, which is needed to calculate the Lake Erie water balance 
C  term, ermwb, below when stage discharge relationship is based on more than
C  one gauge (when erb > 0).
    
      IF (erb .gt. 0.) THEN
         IF (i.eq.1) THEN
            deltaz4(qm,mon) = niaqlev(qm,mon)-niz
         ELSE
            deltaz4(qm,mon) = 0.
         ENDIF
      ENDIF

C  WB's are "water balance" summations...NBS, diversions, consumptive 
C  use, inflow, outflow, and ice retardation. These terms are needed
C  to find the matrix solution in the next subroutine.

      mhqwb = mhqnbs(qm,mon)+miqdiv(qm,mon)+huqdiv(qm,mon)+miqgdw(qm,mon
     &  )+huqgdw(qm,mon)-miqcon(qm,mon)-huqcon(qm,mon)+spqflo(qm,mon)-
     &  outmh+mhqice(qm,mon)

      scqwb = scqnbs(qm,mon)+scqdiv(qm,mon)+scqgdw(qm,mon)-
     &   scqcon(qm,mon)+outmh-mhqice(qm,mon)-outsc+scqice(qm,mon)

      erqwb = erqnbs(qm,mon)+erqdiv(qm,mon)+erqgdw(qm,mon)-
     &   erqcon(qm,mon)+outsc-scqice(qm,mon)-outer+erqice(qm,mon)+
     &   (.5*erb*out6-.5*era*(1.-erwt)*out5)*deltaz4(qm,mon)

      RETURN
      END
C----------------------------------------------------------------------
      Subroutine MATRIX(mhz,scz,erz)  
C     ------------------------------      
C
C  This subroutine uses the coefficients calculated above to solve for the        
C  unknowns deltaz1, deltaz2, and deltaz3 by way of a matrix.  These deltaz's 
C  define the new end-of-increment levels.  This subroutine
C  was written by Anne H. Clites, NOAA/GLERL, November, 1995.
C
      IMPLICIT NONE
      INTEGER qm,mon,inc,nyr
      CHARACTER title(3)*80,ext*4,type*4,units*7
      REAL deltaz4(4,12)
      REAL dt
      REAL a1m1,a2m1,a3m1,b1m1,b2m1,b3m1
      REAL mhz,scz,erz,mhz2,scz2,erz2
      REAL*8 cdet,z1det,z2det,z3det
      REAL*8 deltaz1,deltaz2,deltaz3
      REAL mhqnbs(4,12),scqnbs(4,12),erqnbs(4,12)
      REAL miqpre(4,12),miqrun(4,12),miqevp(4,12),miqcon(4,12)
      REAL huqpre(4,12), huqrun(4,12), huqevp(4,12), huqcon(4,12)
      REAL scqpre(4,12), scqrun(4,12), scqevp(4,12), scqcon(4,12)
      REAL erqpre(4,12), erqrun(4,12), erqevp(4,12), erqcon(4,12)
      REAL miqgdw(4,12), huqgdw(4,12), scqgdw(4,12), erqgdw(4,12)
      REAL miqdiv(4,12), huqdiv(4,12), scqdiv(4,12), erqdiv(4,12)
      REAL mhqice(4,12), scqice(4,12), erqice(4,12), spqflo(4,12)
      REAL mhmnbs(12,100),scmnbs(12,100),ermnbs(12,100)
      REAL*8 c11,c12,c13,c21,c22,c23,c31,c32,c33
      REAL*8 mhqwb,scqwb,erqwb


      COMMON /COEFFS/ c11,c12,c13,c21,c22,c23,c31,c32,c33
      COMMON /EXP/    a1m1,a2m1,a3m1,b1m1,b2m1,b3m1
      COMMON /HYDRO/    miqpre,huqpre,scqpre,erqpre,miqrun,huqrun,
     &                  scqrun,erqrun,miqdiv,huqdiv,scqdiv,erqdiv,
     &                  miqgdw,huqgdw,scqgdw,erqgdw,miqevp,huqevp,
     &                  scqevp,erqevp,miqcon,huqcon,scqcon,erqcon,
     &                  spqflo,mhqice,scqice,erqice,mhqnbs,scqnbs,
     &                  erqnbs,deltaz4,mhmnbs,scmnbs,ermnbs
      COMMON /MISC/    title,ext,type,dt,mon,nyr,inc,qm,units
      COMMON /NEWLEV/ mhz2,scz2,erz2
      COMMON /WATBAL/ mhqwb,scqwb,erqwb
      COMMON /DELTAS/ deltaz1,deltaz2,deltaz3

      mhz2=0.
      scz2=0.
      erz2=0.

C  cdet is the determinant of the C matrix.
C  z1det, z2det, AND z3det are the determinants of the Michigan-Huron,
C  St. Clair, and Erie matrices, respectively.

      cdet = c11*((c22*c33)-(c23*c32)) - c12*((c21*c33)-(c23*c31))
     &       + c13*((c21*c32)-(c22*c31))

      z1det = mhqwb*((c22*c33)-(c23*c32)) - c12*((scqwb*c33)-
     &        (c23*erqwb)) + c13*((scqwb*c32)-(c22*erqwb))

      z2det = c11*((scqwb*c33)-(c23*erqwb)) - mhqwb*((c21*c33)-
     &        (c23*c31)) + c13*((c21*erqwb)-(scqwb*c31))

      z3det = c11*((c22*erqwb)-(scqwb*c32)) - c12*((c21*erqwb)-
     &        (scqwb*c31)) + mhqwb*((c21*c32)-(c22*c31))

C  The change in level equals the determinant divided by the cdet.

      deltaz1 = z1det/cdet     
      deltaz2 = z2det/cdet   
      deltaz3 = z3det/cdet

C  New end-of-increment levels.

      mhz2 = mhz + deltaz1
      scz2 = scz + deltaz2
      erz2 = erz + deltaz3

      RETURN
      END

C------------------------------------------------------------------------------
      Subroutine STATS
C     ----------------
C  This subroutine uses monthly parameters calculated in Main to calculate
C  monthly and annual means, standard deviations, max's and min's. This subroutine
C  was written by Anne Clites and modified by Tim Hunter, June, 1997.
C
      IMPLICIT  NONE
      REAL onminflo(12,100),onqinflo(4,12,100)
      REAL mhmmlev(12,100),mhmmflo(12,100),mhmmbom(12,100)
      REAL scmmlev(12,100),scmmflo(12,100),scmmbom(12,100)
      REAL ermmlev(12,100),ermmflo(12,100),ermmbom(12,100)
      REAL mhlevsum(100),mhbomsum(100),mhflosum(100),sclevsum(100)
      REAL scbomsum(100),scflosum(100),erlevsum(100),erbomsum(100)
      REAL erflosum(100)
      REAL amhlev(12),amhflo(12),amhbom(12),asclev(12),ascflo(12),
     &  ascbom(12),aerlev(12),aerflo(12),aerbom(12)
      REAL zmhminle(12),zmhmaxle(12),zmhminfl(12),zmhmaxfl(12),
     &  zmhminbo(12),zmhmaxbo(12),zscminle(12),zscmaxle(12),
     &  zscminfl(12),zscmaxfl(12),zscminbo(12),zscmaxbo(12),
     &  zerminle(12),zermaxle(12),zerminfl(12),zermaxfl(12),
     &  zerminbo(12),zermaxbo(12)
      REAL mhlevsdv(12),mhbomsdv(12),mhflosdv(12),sclevsdv(12),
     &  scbomsdv(12),scflosdv(12),erlevsdv(12),erbomsdv(12),erflosdv(12)
      REAL aonflo(12),onflosdv(12),zonmaxfl(12),zonminfl(12)
      REAL onflosum(100)
      REAL sdvmhlev, sdvmhbom, sdvmhflo, sdvsclev, sdvscbom, sdvscflo
      REAL sdverlev, sdverbom, sdverflo, sdvonflo
      REAL amaxmhle, amaxmhbo, amaxmhfl, aminmhle, aminmhbo, aminmhfl
      REAL amaxscle, amaxscbo, amaxscfl, aminscle, aminscbo, aminscfl
      REAL amaxerle, amaxerbo, amaxerfl, aminerle, aminerbo, aminerfl
      REAL amaxonfl, aminonfl
      REAL mhlevann, sclevann, erlevann, mhlevmax, sclevmax, erlevmax
      REAL mhlevmin, sclevmin, erlevmin, mhbomann, scbomann, erbomann
      REAL mhbommax, scbommax, erbommax, mhbommin, scbommin, erbommin
      REAL mhfloann, scfloann, erfloann, mhflomax, scflomax, erflomax
      REAL mhflomin, scflomin, erflomin, onfloann, onflomax, onflomin


C  Common block declarations.
                  

      COMMON /LEVFLO/   mhmmlev,scmmlev,ermmlev,mhmmbom,scmmbom,ermmbom,
     &                  mhmmflo,scmmflo,ermmflo,onminflo,onqinflo,
     &                  mhbomsum,mhflosum,sclevsum,scbomsum,scflosum,
     &                  erlevsum,erbomsum,erflosum,mhlevsum,onflosum
      COMMON /STATS1/   amhbom,amhlev,amhflo,ascbom,asclev,ascflo,
     &                  aerbom,aerlev,aerflo,zmhmaxle,zmhminle,zmhmaxbo,
     &                  zmhminbo,zmhmaxfl,zmhminfl,zscmaxle,zscminle,
     &                  zscmaxbo,zscminbo,zscmaxfl,zscminfl,zermaxle,
     &                  zerminle,zermaxbo,zerminbo,zermaxfl,zerminfl
      COMMON /STATS2/   mhlevsdv,mhbomsdv,
     &                  mhflosdv,sclevsdv,scbomsdv,scflosdv,erlevsdv,
     &                  erbomsdv,erflosdv,sdvmhlev,sdvmhbom,sdvmhflo,
     &                  sdvsclev,sdvscbom,sdvscflo,sdverlev,sdverbom,
     &                  sdverflo,sdvonflo,mhlevann,mhbomann,mhfloann,
     &                  scbomann,scfloann,erlevann,erbomann,erfloann,
     &                  mhlevmax,mhlevmin,mhbommax,mhbommin,mhflomax,
     &                  mhflomin,sclevmax,sclevmin,scbommax,scbommin,
     &                  scflomax,scflomin,erlevmax,erlevmin,erbommax,
     &                  erbommin,erflomax,erflomin,aonflo,onfloann,
     &                  onflosdv,zonmaxfl,zonminfl,onflomax,onflomin,
     &                  sclevann
C
C  Compute statistics for Michigan-Huron mean levels, flows, beginning-of-month levels.
C
      CALL STATSUB(mhlevsum, mhmmlev, mhlevann, mhlevmin, mhlevmax,
     &             sdvmhlev, amhlev,  zmhminle, zmhmaxle, aminmhle,
     &             amaxmhle, mhlevsdv)

      CALL STATSUB(mhflosum, mhmmflo, mhfloann, mhflomin, mhflomax,
     &             sdvmhflo, amhflo,  zmhminfl, zmhmaxfl, aminmhfl,
     &             amaxmhfl, mhflosdv)

      CALL STATSUB(mhbomsum, mhmmbom, mhbomann, mhbommin, mhbommax,
     &             sdvmhbom, amhbom,  zmhminbo, zmhmaxbo, aminmhbo,
     &             amaxmhbo, mhbomsdv)

C
C  Compute statistics for St. Clair mean levels, flows, beginning-of-month levels.
C
      CALL STATSUB(sclevsum, scmmlev, sclevann, sclevmin, sclevmax,
     &             sdvsclev, asclev,  zscminle, zscmaxle, aminscle,
     &             amaxscle, sclevsdv)

      CALL STATSUB(scflosum, scmmflo, scfloann, scflomin, scflomax,
     &             sdvscflo, ascflo,  zscminfl, zscmaxfl, aminscfl,
     &             amaxscfl, scflosdv)

      CALL STATSUB(scbomsum, scmmbom, scbomann, scbommin, scbommax,
     &             sdvscbom, ascbom,  zscminbo, zscmaxbo, aminscbo,
     &             amaxscbo, scbomsdv)

C
C  Compute statistics for Erie mean levels, flows, beginning-of-month levels.
C
      CALL STATSUB(erlevsum, ermmlev, erlevann, erlevmin, erlevmax,
     &             sdverlev, aerlev,  zerminle, zermaxle, aminerle,
     &             amaxerle, erlevsdv)

      CALL STATSUB(erflosum, ermmflo, erfloann, erflomin, erflomax,
     &             sdverflo, aerflo,  zerminfl, zermaxfl, aminerfl,
     &             amaxerfl, erflosdv)

      CALL STATSUB(erbomsum, ermmbom, erbomann, erbommin, erbommax,
     &             sdverbom, aerbom,  zerminbo, zermaxbo, aminerbo,
     &             amaxerbo, erbomsdv)

C
C  Compute statistics for Ontario flows
C
      CALL STATSUB(onflosum, onminflo, onfloann, onflomin, onflomax,
     &             sdvonflo, aonflo,  zonminfl, zonmaxfl, aminonfl,
     &             amaxonfl, onflosdv)

      RETURN
      END

C----------------------------------------------------------------------
      SUBROUTINE STATSUB(annmeanin,   monmeanin, 
     &                   annualmean,  annualmin,  annualmax,  annsd,
     &                   monthlymean, monthlymin, monthlymax,
     &                   monavgmin,   monavgmax,  monsd)

C
C  This subroutine was written by Tim Hunter in June, 1997.
C

      IMPLICIT  NONE
      REAL      annmeanin(100), monmeanin(12,100), annsd
      REAL      annualmean, annualmin, annualmax
      REAL      monthlymean(12), monthlymin(12), monthlymax(12)
      REAL      monavgmin, monavgmax, monsd(12)
      REAL      total, mtotal(12)

      INTEGER   iyr1, mon1, iyrn, monn
      INTEGER   m, y, numyears, mn1, mn2, mcount(12)

      COMMON /PERIODS/ iyr1, mon1, iyrn, monn

C
C  Compute "long-term" annual mean, min, and max values.
C
      numyears = iyrn - iyr1 + 1 
      annualmin = annmeanin(1)
      annualmax = annmeanin(1)
      total = 0.
      DO y = 1, numyears
         total = total + annmeanin(y)
         annualmin = MIN(annualmin, annmeanin(y))
         annualmax = MAX(annualmax, annmeanin(y))
      ENDDO
      annualmean = total / numyears

C
C  Compute the standard deviation of the series of annual means:
C     1) sum the squared differences of each annual mean from the long-term mean
C     2) Std Dev = the square root of the variance.
C
      total = 0.
      DO y = 1, numyears
         total = total + (annmeanin(y) - annualmean)**2
      ENDDO
      annsd = SQRT(total / numyears)

C
C  Calculate "long-term" monthly mean, min, max values for each month
C      (e.g. min January value, max Jan value, mean Jan value)
C
      DO m = 1, 12
        mtotal(m) = 0.
        mcount(m) = 0
      ENDDO

C
C  Initialize the max and min arrays.
C  If there is only one year of interest, then I
C  only want to assign values for months 'mon1' to 'monn'.  If there are
C  2 or more years, then we will do it for all 12 months of the year.
C  The part that makes it messy is that we need to make sure the value we
C  assign is valid; that is, that it is one of the months for which we
C  have legitimate data. So we can't just blindly do months 1 to 12.
C  If we did, then the following cases apply:
C    1) Only one year
C         Months 1 to mon1-1 would be assigned zero.
C         Months monn+1 to 12 would be assigned zero.
C         This would mess up the minimum values, for sure.
C    2) Only 2 years
C         Can't just use months 1-12 of year 2 as initial values
C         because then months monn+1 to 12 would be assigned zero.
C    3) More than 2 years
C         We could legitimately use months 1-12 of year 2, but I think
C         it is "cleaner" not to search for this special case.
C
C  Since we want to handle the general cases as well as the special case,
C  we will simply ignore years 3 to n.  Initial values will come from
C  years 1 and/or 2 only.
C
      DO y = 1, numyears
         IF (y .LE. 2) THEN
           mn1 = 1
           mn2 = 12
           IF (y .EQ. 1)        mn1 = mon1
           IF (y .EQ. numyears) mn2 = monn
           DO m = mn1, mn2
             monthlymin(m) = monmeanin(m,y)
             monthlymax(m) = monmeanin(m,y)
           ENDDO
         ENDIF
      ENDDO

      DO y = 1, numyears
         mn1 = 1
         mn2 = 12
         IF (y .EQ. 1)        mn1 = mon1
         IF (y .EQ. numyears) mn2 = monn
         DO m = mn1, mn2
            mtotal(m) = mtotal(m) + monmeanin(m,y)
            mcount(m) = mcount(m) + 1
            monthlymin(m) = MIN(monthlymin(m), monmeanin(m,y))
            monthlymax(m) = MAX(monthlymax(m), monmeanin(m,y))
         ENDDO
      ENDDO

      DO m = 1, 12
         monthlymean(m) = -999999.9
         IF (mcount(m) .GT. 0) monthlymean(m) = mtotal(m) / mcount(m)
      ENDDO

C
C  Compute the standard deviation of each series of monthly means
C
C     1) Sum the squared differences of each monthly mean from the long-term mean
C        for that month of the year
C     2) Std Dev = the square root of the variance
C
      DO m = 1, 12
         mtotal(m) = 0.
         mcount(m) = 0
      ENDDO

      DO y = 1, numyears
         mn1 = 1
         mn2 = 12
         IF (y .EQ. 1)        mn1 = mon1
         IF (y .EQ. numyears) mn2 = monn
         DO m = mn1, mn2
            mtotal(m) = mtotal(m) + (monmeanin(m,y)-monthlymean(m))**2
            mcount(m) = mcount(m) + 1
         ENDDO
      ENDDO

      DO m = 1, 12
         monsd(m) = -999999.9
         IF (mcount(m) .GT. 0) monsd(m) = SQRT(mtotal(m) / mcount(m))
      ENDDO

C
C  Determine the max and min monthly average value.
C  i.e., which of the 12 monthly average values is the highest and lowest?
C
      mn1 = 1
      mn2 = 12
      IF (numyears .EQ. 1) THEN
         mn1 = mon1
         mn2 = monn
      ENDIF
      monavgmax = monthlymean(mn1)
      monavgmin = monthlymean(mn1)
      DO m = 1, 12
         IF (monthlymean(m) .GT. -999999.8) THEN
            monavgmax = MAX(monavgmax, monthlymean(m))
            monavgmin = MIN(monavgmin, monthlymean(m))
         ENDIF
      ENDDO

      RETURN
      END

C------------------------------------------------------------------------------
      Subroutine OUTDATA
C     -------------------
C      
C  Outdata opens output files and prints to them. This subroutine was written by
C  Anne H. Clites in February, 1997. 
      
      IMPLICIT NONE
      CHARACTER fnames(34)*40,title(3)*80,ext*4,type*4,months(12)*9,
     &  dat*8,tim*8,stat(4)*4,run*3,chyear*4,units*7
      INTEGER  numyears, totmonths, m1, m2, m3, m4, nfile, nruns
      INTEGER  iyr1,iyrn,mon1,monn,k,nyr,nyear,i,iq,inc,mon,qm
      REAL mhk,mhym,mha,mhb,mhwt,dt,mhsar,misar,husar
      REAL sck,scym,sca,scb,scwt,scsar
      REAL erk,erym,era,erb,erwt,ersar
      REAL mhmmlev(12,100),mhmmflo(12,100),mhmmbom(12,100)
      REAL scmmlev(12,100),scmmflo(12,100),scmmbom(12,100)
      REAL ermmlev(12,100),ermmflo(12,100),ermmbom(12,100)
      REAL mhlevsum(100),mhbomsum(100),mhflosum(100),sclevsum(100)
      REAL scbomsum(100),scflosum(100),erlevsum(100),erbomsum(100)
      REAL erflosum(100)
      REAL mhlevmax,mhlevmin,mhlevann,mhbommax,mhbommin,mhbomann
      REAL mhflomax,mhflomin,mhfloann
      REAL sdvmhlev,sdvmhbom,sdvmhflo,sdvsclev,sdvscbom,sdvscflo
      REAL sdverlev,sdverbom,sdverflo
      REAL sclevann,scbomann,scfloann,sclevmax,sclevmin,scbommax
      REAL scbommin,scflomax,scflomin,onfloann,onflomax,onflomin
      REAL sdvonflo,erlevann,erbomann,erfloann,erlevmax,erlevmin
      REAL erbommax,erflomax,erflomin,erbommin      
      REAL amhlev(12),amhflo(12),amhbom(12),asclev(12),ascflo(12),
     &  ascbom(12),aerlev(12),aerflo(12),aerbom(12)
      REAL zmhminle(12),zmhmaxle(12),zmhminfl(12),zmhmaxfl(12),
     &  zmhminbo(12),zmhmaxbo(12),zscminle(12),zscmaxle(12),
     &  zscminfl(12),zscmaxfl(12),zscminbo(12),zscmaxbo(12),
     &  zerminle(12),zermaxle(12),zerminfl(12),zermaxfl(12),
     &  zerminbo(12),zermaxbo(12)
      REAL mhlevsdv(12),mhbomsdv(12),mhflosdv(12),sclevsdv(12),
     &  scbomsdv(12),scflosdv(12),erlevsdv(12),erbomsdv(12),erflosdv(12)
      REAL aonflo(12),onflosdv(12),zonmaxfl(12),zonminfl(12)
      REAL onflosum(100)
      REAL miqpre(4,12), miqrun(4,12), miqevp(4,12), miqcon(4,12)
      REAL huqpre(4,12), huqrun(4,12), huqevp(4,12), huqcon(4,12)
      REAL scqpre(4,12), scqrun(4,12), scqevp(4,12), scqcon(4,12)
      REAL erqpre(4,12), erqrun(4,12), erqevp(4,12), erqcon(4,12)
      REAL miqgdw(4,12), huqgdw(4,12), scqgdw(4,12), erqgdw(4,12)
      REAL miqdiv(4,12), huqdiv(4,12), scqdiv(4,12), erqdiv(4,12)
      REAL mhqice(4,12), scqice(4,12), erqice(4,12), spqflo(4,12)
      REAL mhmnbs(12,100),scmnbs(12,100),ermnbs(12,100)
      REAL mhqnbs(4,12),scqnbs(4,12),erqnbs(4,12)
      REAL mhz,scz,erz,niz,deltaz4(4,12)
      REAL onminflo(12,100),onqinflo(4,12,100)

C  Common block declarations.
                  
      COMMON /BOP/      mhz,scz,erz,niz
      COMMON /FILENAME/ fnames,nfile,nruns,run
      COMMON /HYDRO/    miqpre,huqpre,scqpre,erqpre,miqrun,huqrun,
     &                  scqrun,erqrun,miqdiv,huqdiv,scqdiv,erqdiv,
     &                  miqgdw,huqgdw,scqgdw,erqgdw,miqevp,huqevp,
     &                  scqevp,erqevp,miqcon,huqcon,scqcon,erqcon,
     &                  spqflo,mhqice,scqice,erqice,mhqnbs,scqnbs,
     &                  erqnbs,deltaz4,mhmnbs,scmnbs,ermnbs
      COMMON /LEVFLO/   mhmmlev,scmmlev,ermmlev,mhmmbom,scmmbom,ermmbom,
     &                  mhmmflo,scmmflo,ermmflo,onminflo,onqinflo,
     &                  mhbomsum,mhflosum,sclevsum,scbomsum,scflosum,
     &                  erlevsum,erbomsum,erflosum,mhlevsum,onflosum
      COMMON /MISC/     title,ext,type,dt,mon,nyr,inc,qm,units
      COMMON /PERIODS/  iyr1,mon1,iyrn,monn
      COMMON /SA/       mhsar,misar,husar,scsar,ersar
      COMMON /STAGE/    mhk,sck,erk,mhym,scym,erym,mha,sca,era,
     &                  mhb,scb,erb,mhwt,scwt,erwt
      COMMON /STATS1/   amhbom,amhlev,amhflo,ascbom,asclev,ascflo,
     &                  aerbom,aerlev,aerflo,zmhmaxle,zmhminle,zmhmaxbo,
     &                  zmhminbo,zmhmaxfl,zmhminfl,zscmaxle,zscminle,
     &                  zscmaxbo,zscminbo,zscmaxfl,zscminfl,zermaxle,
     &                  zerminle,zermaxbo,zerminbo,zermaxfl,zerminfl
      COMMON /STATS2/   mhlevsdv,mhbomsdv,
     &                  mhflosdv,sclevsdv,scbomsdv,scflosdv,erlevsdv,
     &                  erbomsdv,erflosdv,sdvmhlev,sdvmhbom,sdvmhflo,
     &                  sdvsclev,sdvscbom,sdvscflo,sdverlev,sdverbom,
     &                  sdverflo,sdvonflo,mhlevann,mhbomann,mhfloann,
     &                  scbomann,scfloann,erlevann,erbomann,erfloann,
     &                  mhlevmax,mhlevmin,mhbommax,mhbommin,mhflomax,
     &                  mhflomin,sclevmax,sclevmin,scbommax,scbommin,
     &                  scflomax,scflomin,erlevmax,erlevmin,erbommax,
     &                  erbommin,erflomax,erflomin,aonflo,onfloann,
     &                  onflosdv,zonmaxfl,zonminfl,onflomax,onflomin,
     &                  sclevann

      DATA months /'January  ','February ','March    ','April    ',
     &             'May      ','June     ','July     ','August   ',
     &             'September','October  ','November ','December '/
      DATA stat  /'Mean','StDv','Max ','Min '/


C
C  If some of the statistical values were never given valid values,
C  then assign the missing data code to them (-10000000.)
C
      totmonths = (iyrn-iyr1)*12 + monn-mon1
      IF (totmonths .LT. 12) THEN
         m1 = 0
         m2 = -1
         m3 = 0
         m4 = -1
         IF (mon1 .GT. monn) THEN
            m1 = monn + 1
            m2 = mon1 - 1
         ELSE
            IF (mon1 .GT. 1) THEN
               m1 = 1
               m2 = mon1 - 1
            ENDIF
            IF (monn .LT. 12) THEN
               m3 = monn + 1
               m4 = 12
            ENDIF
         ENDIF
         DO k = m1, m2
            amhbom(k) = -10000000.
            amhlev(k) = -10000000.
            amhflo(k) = -10000000.
            ascbom(k) = -10000000.
            asclev(k) = -10000000.
            ascflo(k) = -10000000.
            aerbom(k) = -10000000.
            aerlev(k) = -10000000.
            aerflo(k) = -10000000.
            mhbomsdv(k) = -10000000.
            mhlevsdv(k) = -10000000.
            mhflosdv(k) = -10000000.
            scbomsdv(k) = -10000000.
            sclevsdv(k) = -10000000.
            scflosdv(k) = -10000000.
            erbomsdv(k) = -10000000.
            erlevsdv(k) = -10000000.
            erflosdv(k) = -10000000.
            zmhmaxbo(k) = -10000000.
            zmhmaxle(k) = -10000000.
            zmhmaxfl(k) = -10000000.
            zscmaxbo(k) = -10000000.
            zscmaxle(k) = -10000000.
            zscmaxfl(k) = -10000000.
            zermaxbo(k) = -10000000.
            zermaxle(k) = -10000000.
            zermaxfl(k) = -10000000.
            zmhminbo(k) = -10000000.
            zmhminle(k) = -10000000.
            zmhminfl(k) = -10000000.
            zscminbo(k) = -10000000.
            zscminle(k) = -10000000.
            zscminfl(k) = -10000000.
            zerminbo(k) = -10000000.
            zerminle(k) = -10000000.
            zerminfl(k) = -10000000.
         ENDDO
         DO k = m3, m4
            amhbom(k) = -10000000.
            amhlev(k) = -10000000.
            amhflo(k) = -10000000.
            ascbom(k) = -10000000.
            asclev(k) = -10000000.
            ascflo(k) = -10000000.
            aerbom(k) = -10000000.
            aerlev(k) = -10000000.
            aerflo(k) = -10000000.
            mhbomsdv(k) = -10000000.
            mhlevsdv(k) = -10000000.
            mhflosdv(k) = -10000000.
            scbomsdv(k) = -10000000.
            sclevsdv(k) = -10000000.
            scflosdv(k) = -10000000.
            erbomsdv(k) = -10000000.
            erlevsdv(k) = -10000000.
            erflosdv(k) = -10000000.
            zmhmaxbo(k) = -10000000.
            zmhmaxle(k) = -10000000.
            zmhmaxfl(k) = -10000000.
            zscmaxbo(k) = -10000000.
            zscmaxle(k) = -10000000.
            zscmaxfl(k) = -10000000.
            zermaxbo(k) = -10000000.
            zermaxle(k) = -10000000.
            zermaxfl(k) = -10000000.
            zmhminbo(k) = -10000000.
            zmhminle(k) = -10000000.
            zmhminfl(k) = -10000000.
            zscminbo(k) = -10000000.
            zscminle(k) = -10000000.
            zscminfl(k) = -10000000.
            zerminbo(k) = -10000000.
            zerminle(k) = -10000000.
            zerminfl(k) = -10000000.
         ENDDO
      ENDIF

      numyears = iyrn - iyr1 + 1
      DO nyr = 1, numyears
         nyear = iyr1 + nyr - 1
         m1 = 1
         m2 = 12
         IF (nyr .EQ. 1) THEN
            m2 = mon1 - 1
            DO i = m1, m2
               DO iq = 1, 4
                 onqinflo(iq,i,nyr) = -10000000.
               ENDDO
            ENDDO
         ENDIF
         IF (nyr .EQ. numyears) THEN
            m1 = monn + 1
            DO i = m1, m2
               DO iq = 1, 4
                 onqinflo(iq,i,nyr) = -10000000.
               ENDDO
            ENDDO
         ENDIF
      ENDDO

C  Open output files. The 'summary' file is a look at levels and
C  flows.  For each lake, there is a separate file for mean levels,
C  beginning-of-month levels, and mean flows. The 'details' file is a
C  more lengthy summary file containing all inputs. All output files include
C  some brief statistics (AVG, STD, MAX, MIN) for monthly and annual values.
                        
      OPEN(UNIT=8, FILE='..\outdata\summary'//ext, STATUS='UNKNOWN')
      OPEN(UNIT=9, FILE='..\outdata\onqinflo'//ext, STATUS='UNKNOWN')
      OPEN(UNIT=11, FILE='..\outdata\mhmlev'//ext, STATUS='UNKNOWN')
      OPEN(UNIT=12, FILE='..\outdata\mhblev'//ext, STATUS='UNKNOWN')
      OPEN(UNIT=13, FILE='..\outdata\mhmflo'//ext, STATUS='UNKNOWN')
      OPEN(UNIT=21, FILE='..\outdata\scmlev'//ext, STATUS='UNKNOWN')
      OPEN(UNIT=22, FILE='..\outdata\scblev'//ext, STATUS='UNKNOWN')
      OPEN(UNIT=23, FILE='..\outdata\scmflo'//ext, STATUS='UNKNOWN')
      OPEN(UNIT=31, FILE='..\outdata\ermlev'//ext, STATUS='UNKNOWN')
      OPEN(UNIT=32, FILE='..\outdata\erblev'//ext, STATUS='UNKNOWN')
      OPEN(UNIT=33, FILE='..\outdata\ermflo'//ext, STATUS='UNKNOWN')
      OPEN(UNIT=34, FILE='..\outdata\onminflo'//ext, STATUS='UNKNOWN')

C  Write header information to output files.
      
      CALL DATE(dat)
      CALL TIME(tim)
      WRITE(8,1013) 'summary'//ext,dat,tim
      WRITE(8,1000)
      WRITE(8,1001) (title(i),i=1,3)
      WRITE(8,1002) inc
      WRITE(8,1021) units
      WRITE(8,1003) mhk,mhym,mha,mhb,mhwt
      WRITE(8,1004) sck,scym,sca,scb,scwt
      WRITE(8,1005) erk,erym,era,erb,erwt
      WRITE(8,1007)
      WRITE(9,1013) 'onqinflo'//ext,dat,tim
      WRITE(9,1200)
      WRITE(11,1013) 'mhmlev'//ext,dat,tim
      WRITE(11,2011)
      WRITE(12,1013) 'mhblev'//ext,dat,tim
      WRITE(12,2012)
      WRITE(13,1013) 'mhmflo'//ext,dat,tim
      WRITE(13,2013)
      WRITE(21,1013) 'scmlev'//ext,dat,tim
      WRITE(21,2021)
      WRITE(22,1013) 'scblev'//ext,dat,tim
      WRITE(22,2022)
      WRITE(23,1013) 'scmflo'//ext,dat,tim
      WRITE(23,2023)
      WRITE(31,1013) 'ermlev'//ext,dat,tim
      WRITE(31,2031)
      WRITE(32,1013) 'erblev'//ext,dat,tim
      WRITE(32,2032)
      WRITE(33,1013) 'ermflo'//ext,dat,tim
      WRITE(33,2033)
      WRITE(34,1013) 'onminflo'//ext,dat,tim
      WRITE(34,1202)

C  Write monthly results to both summary file and individual files.

      numyears = iyrn - iyr1 + 1
      DO nyr = 1, numyears
         nyear = iyr1 + nyr - 1
         m1 = 1
         m2 = 12
         IF (nyr .EQ. 1)        m1 = mon1
         IF (nyr .EQ. numyears) m2 = monn

C   Print monthly values to summary file...one line per month.

         DO mon = m1, m2
            WRITE(8,1010) nyear,mon,mhmnbs(mon,nyr),mhmmbom(mon,nyr),
     &              mhmmlev(mon,nyr),mhmmflo(mon,nyr),scmnbs(mon,nyr),
     &              scmmbom(mon,nyr),scmmlev(mon,nyr),scmmflo(mon,nyr),
     &              ermnbs(mon,nyr),ermmbom(mon,nyr),ermmlev(mon,nyr), 
     &              ermmflo(mon,nyr)
         ENDDO

C  Write out quarter-monthly Ontario inflows (Erie outflows + Welland
C  diversion).  Also write it out monthly for plotting
C  to compare with BOC plots.

         DO iq=1,4
              WRITE(9,1201) nyear,iq,(onqinflo(iq,mon,nyr),mon=1,12)
         ENDDO

C  Print monthlies to individual files...one line per year.
  
         WRITE(chyear, '(I4)') nyear

         CALL MONOUTLINE(11,2,chyear,mhmmlev,mhlevsum,nyr,m1,m2)
         CALL MONOUTLINE(12,2,chyear,mhmmbom,mhbomsum,nyr,m1,m2)
         CALL MONOUTLINE(13,0,chyear,mhmmflo,mhflosum,nyr,m1,m2)

         CALL MONOUTLINE(21,2,chyear,scmmlev,sclevsum,nyr,m1,m2)
         CALL MONOUTLINE(22,2,chyear,scmmbom,scbomsum,nyr,m1,m2)
         CALL MONOUTLINE(23,0,chyear,scmmflo,scflosum,nyr,m1,m2)

         CALL MONOUTLINE(31,2,chyear,ermmlev,erlevsum,nyr,m1,m2)
         CALL MONOUTLINE(32,2,chyear,ermmbom,erbomsum,nyr,m1,m2)
         CALL MONOUTLINE(33,0,chyear,ermmflo,erflosum,nyr,m1,m2)

         CALL MONOUTLINE(34,0,chyear,onminflo,onflosum,nyr,m1,m2)

         WRITE(8,'(/)')     !  Advance to new page in summary file.
      ENDDO


C  Skip a line in individual files.

      WRITE(9,'(/)')
      WRITE(11,'(/)')
      WRITE(12,'(/)')
      WRITE(13,'(/)')
      WRITE(21,'(/)')
      WRITE(22,'(/)')
      WRITE(23,'(/)')
      WRITE(31,'(/)')
      WRITE(32,'(/)')
      WRITE(33,'(/)')
      WRITE(34,'(/)')

C  Print means for each month and parameter for the period of run.

      WRITE(8,1014)
      WRITE(8,1016)
      
C
C  Set values for the number of years and which months to start/end
C  for statistical outputs
C
      numyears = iyrn - iyr1 + 1
      m1 = 1
      m2 = 12
      IF (numyears .EQ. 1) THEN
         m1 = mon1
         m2 = monn
      ENDIF

C  This prints out the monthly means matrix.

      DO k=m1,m2
         CALL MONSTATSLINE(8, months(k),amhbom(k),amhlev(k),
     &                     amhflo(k),ascbom(k),asclev(k),ascflo(k),
     &                     aerbom(k),aerlev(k),aerflo(k))
      ENDDO

C  Print out the monthly standard deviations matrix.

      WRITE(8,1017)
      WRITE(8,1016)

      DO k=m1,m2
         CALL MONSTATSLINE(8, months(k),mhbomsdv(k),mhlevsdv(k),
     &                     mhflosdv(k),scbomsdv(k),sclevsdv(k),
     &                     scflosdv(k),erbomsdv(k),erlevsdv(k),
     &                     erflosdv(k))
      ENDDO
      

C  Print max and min matrices for monthly parameters for period of run.
  
      WRITE(8,1019)
      WRITE(8,1016)
          
      DO k=m1,m2
         CALL MONSTATSLINE(8, months(k),zmhmaxbo(k),zmhmaxle(k),
     &               zmhmaxfl(k),zscmaxbo(k),zscmaxle(k),zscmaxfl(k),
     &               zermaxbo(k),zermaxle(k),zermaxfl(k))
      ENDDO

      WRITE(8,1020)
      WRITE(8,1016)

      DO k=m1,m2
         CALL MONSTATSLINE(8, months(k),zmhminbo(k),zmhminle(k),
     &               zmhminfl(k),zscminbo(k),zscminle(k),zscminfl(k),
     &               zerminbo(k),zerminle(k),zerminfl(k))
      ENDDO


C  Print mean, standard deviation, max, and min to individual files.

      CALL STATOUTLINE(11, 2, stat(1), amhlev,   mhlevann)
      CALL STATOUTLINE(11, 2, stat(2), mhlevsdv, sdvmhlev)
      CALL STATOUTLINE(11, 2, stat(3), zmhmaxle, mhlevmax)
      CALL STATOUTLINE(11, 2, stat(4), zmhminle, mhlevmin)

      CALL STATOUTLINE(12, 2, stat(1), amhbom,   mhbomann)
      CALL STATOUTLINE(12, 2, stat(2), mhbomsdv, sdvmhbom)
      CALL STATOUTLINE(12, 2, stat(3), zmhmaxbo, mhbommax)
      CALL STATOUTLINE(12, 2, stat(4), zmhminbo, mhbommin)
      
      CALL STATOUTLINE(13, 0, stat(1), amhflo,   mhfloann)
      CALL STATOUTLINE(13, 0, stat(2), mhflosdv, sdvmhflo)
      CALL STATOUTLINE(13, 0, stat(3), zmhmaxfl, mhflomax)
      CALL STATOUTLINE(13, 0, stat(4), zmhminfl, mhflomin)
      
      CALL STATOUTLINE(21, 2, stat(1), asclev,   sclevann)
      CALL STATOUTLINE(21, 2, stat(2), sclevsdv, sdvsclev)
      CALL STATOUTLINE(21, 2, stat(3), zscmaxle, sclevmax)
      CALL STATOUTLINE(21, 2, stat(4), zscminle, sclevmin)

      CALL STATOUTLINE(22, 2, stat(1), ascbom,   scbomann)
      CALL STATOUTLINE(22, 2, stat(2), scbomsdv, sdvscbom)
      CALL STATOUTLINE(22, 2, stat(3), zscmaxbo, scbommax)
      CALL STATOUTLINE(22, 2, stat(4), zscminbo, scbommin)
      
      CALL STATOUTLINE(23, 0, stat(1), ascflo,   scfloann)
      CALL STATOUTLINE(23, 0, stat(2), scflosdv, sdvscflo)
      CALL STATOUTLINE(23, 0, stat(3), zscmaxfl, scflomax)
      CALL STATOUTLINE(23, 0, stat(4), zscminfl, scflomin)
      
      CALL STATOUTLINE(31, 2, stat(1), aerlev,   erlevann)
      CALL STATOUTLINE(31, 2, stat(2), erlevsdv, sdverlev)
      CALL STATOUTLINE(31, 2, stat(3), zermaxle, erlevmax)
      CALL STATOUTLINE(31, 2, stat(4), zerminle, erlevmin)

      CALL STATOUTLINE(32, 2, stat(1), aerbom,   erbomann)
      CALL STATOUTLINE(32, 2, stat(2), erbomsdv, sdverbom)
      CALL STATOUTLINE(32, 2, stat(3), zermaxbo, erbommax)
      CALL STATOUTLINE(32, 2, stat(4), zerminbo, erbommin)
      
      CALL STATOUTLINE(33, 0, stat(1), aerflo,   erfloann)
      CALL STATOUTLINE(33, 0, stat(2), erflosdv, sdverflo)
      CALL STATOUTLINE(33, 0, stat(3), zermaxfl, erflomax)
      CALL STATOUTLINE(33, 0, stat(4), zerminfl, erflomin)
      
      CALL STATOUTLINE(34, 0, stat(1), aonflo,   onfloann)
      CALL STATOUTLINE(34, 0, stat(2), onflosdv, sdvonflo)
      CALL STATOUTLINE(34, 0, stat(3), zonmaxfl, onflomax)
      CALL STATOUTLINE(34, 0, stat(4), zonminfl, onflomin)
      

C  Format statements.

 1000 FORMAT (/15X,'MIDLAKES - Middle Lakes Hydrologic Response Model'/)   
 1001 FORMAT ( A80 )   
 1002 FORMAT (/'Time increments per month:',I3/ )            
 1003 FORMAT ('Stage-discharge equation parameters:'/20X,'-K-',7X,'-YM-'
     &,5X,'-A-',4X,'-B-',4X,'-WT-',
     &/,'Michigan-Huron: ',F9.4,4X,F6.2,4X,F3.1,4X,F3.1,4X,F3.1)
 1004 FORMAT ('St. Clair: '5X,F9.4,4X,F6.2,4X,F3.1,4X,F3.1,4X,F3.1)
 1005 FORMAT ('Erie: ',10X,F9.4,4X,F6.2,4X,F3.1,4X,F3.1,4X,F3.1,/)
 1006 FORMAT (I4,3X,3(1X,F7.2,F7.2,1X,F7.0)) 
 1007 FORMAT (10X,'LAKE MICH-HURON',10X,'LAKE ST. CLAIR',12X,'LAKE ERIE'
     &/'Year',1X,'mo.',2X,'NBS',3X,'BOM',3X,'MEAN',4X,'QSC',4X,'NBS',3X,
     &'BOM',3X,'MEAN',4X,'QD',6X,'NBS',3X,'BOM',3X,'MEAN',4X,'QN'/)      
 1008 FORMAT (//33X,'******************'/35X'Annual Means'/33X,
     &'******************'///10X,
     &'LAKE MICH-HURON',10X,'LAKE ST. CLAIR',12X,'LAKE ERIE'/
     &'Year',6X,'BOM',3X,'MEAN',4X,'QSC',6X,'BOM',       
     &3X,'MEAN',4X,'QD',7X,'BOM',3X,'MEAN',4X,'QN'/) 
 1009 FORMAT (//'*** Statistics for Annual Means ***'//10X,'LAKE MICH-
     &HURON',10X,'LAKE ST. CLAIR',12X,'LAKE ERIE'/11X,'BOM',4X,'MEAN',
     &4X,'QSC',6X,'BOM',3X,'MEAN',4X,'QD',7X,'BOM',3X,'MEAN',4X,'QN'/) 
 1010 FORMAT (I4,1X,I2,3(F6.0,F7.2,F7.2,1X,F7.0))  
 1011 FORMAT (7X,3(1X,F7.2,F7.2,1X,F7.0))
 1012 FORMAT (//'*** Statistics for Monthly Means ***'//10x,'LAKE MICH-
     &HURON',10X,'LAKE ST. CLAIR',12X,'LAKE ERIE'/13X,'BOM',4X,'MEAN',
     &4X,'QSC',6X,'BOM',3X,'MEAN',4X,'QD',7X,'BOM',3X,'MEAN',4X,'QN'/) 
 1013 FORMAT (A12,40X,A8,2X,A8)
 1014 FORMAT ('+'//33X,'******************'/35X,'Monthly Means'/
     &33X,'******************'//)
 1015 FORMAT (A9,3(F7.2,1X,F7.2,1X,F7.0))
 1016 FORMAT (10X,'LAKE MICH-HURON',10X,'LAKE ST. CLAIR',12X,'LAKE ERIE'
     &/3X,'Month',3X,'BOM',3X,'MEAN',4X,'QSC',6X,'BOM',
     &3X,'MEAN',4X,'QD',7X,'BOM',3X,'MEAN',4X,'QN'/)      
 1017 FORMAT (//25X,'***************************'/25X,'Monthly Standard 
     &Deviations'/25X,'***************************'//)
 1018 FORMAT (/,A9,3(F7.2,1X,F7.2,1X,F7.0)) 
 1019 FORMAT (//33X,'******************'/34X,'Monthly Maximums'/
     &33X,'******************'//)
 1020 FORMAT (//33X,'******************'/34X,'Monthly Minimums'/
     &33X,'******************'//)
 1021 FORMAT ('Units: ',A7,/)
 1200 FORMAT ('Quarter Monthly Lake Ontario Inflows - (Erie Outflow + We
     &lland Canal)',/,'YEAR',1X,'Q.',2X,'JAN',5X,'FEB',5X,
     &'MAR',5X,'APR',5X,'MAY',5X,'JUN',5X,'JUL',5X,'AUG',5X,'SEP',5X
     &,'OCT',5X,'NOV',5X,'DEC'/)  
 1201 FORMAT (I4,I2,12(F8.0))
 1202 FORMAT ('Monthly Lake Ontario Inflows - (Erie Outflow + Welland Ca
     &nal)',/,'YEAR',3X,'JAN',5X,'FEB',5X,'MAR',5X,'APR',5X,'MAY',5X,
     &'JUN',5X,'JUL',5X,'AUG',5X,'SEP',5X,'OCT',5X,'NOV',5X,'DEC',5X,
     &'ANN',/)  
 1300 FORMAT (3X,I2,F7.3,3X,F7.3,3X,F7.3)                                             
 
 2000 FORMAT (I4, 12F8.2, F8.2) 
 2001 FORMAT (I4, 12F8.0, F8.0)
 2002 FORMAT (A4, 12F8.2, F8.2)
 2003 FORMAT (A4, 12F8.0, F8.0)
 2011 FORMAT (2X,'Lakes Michigan-Huron Monthly Mean Lake Levels',//
     &,'YEAR',3X,'JAN',5X,'FEB',5X,'MAR',5X,'APR',5X,'MAY',5X,'JUN',  
     &5X,'JUL',5X,'AUG',5X,'SEP',5X,'OCT',5X,'NOV',5X,'DEC',5X,'ANN'/)  
 2012 FORMAT (2X,'Lakes Michigan-Huron Beginning-of-Month Lake Levels',/
     &/,'YEAR',3X,'JAN',5X,'FEB',5X,'MAR',5X,'APR',5X,'MAY',5X,'JUN',  
     &5X,'JUL',5X,'AUG',5X,'SEP',5X,'OCT',5X,'NOV',5X,'DEC',5X,'ANN'/)  
 2013 FORMAT (2X,'Lakes Michigan-Huron Monthly Mean Outflows',//
     &,'YEAR',3X,'JAN',5X,'FEB',5X,'MAR',5X,'APR',5X,'MAY',5X,'JUN',  
     &5X,'JUL',5X,'AUG',5X,'SEP',5X,'OCT',5X,'NOV',5X,'DEC',5X,'ANN'/)  
 2021 FORMAT (2X,'Lake St. Clair Monthly Mean Lake Levels',//
     &,'YEAR',3X,'JAN',5X,'FEB',5X,'MAR',5X,'APR',5X,'MAY',5X,'JUN',  
     &5X,'JUL',5X,'AUG',5X,'SEP',5X,'OCT',5X,'NOV',5X,'DEC',5X,'ANN'/)  
 2022 FORMAT (2X,'Lake St. Clair Beginning-of-Month Lake Levels',//
     &,'YEAR',3X,'JAN',5X,'FEB',5X,'MAR',5X,'APR',5X,'MAY',5X,'JUN',  
     &5X,'JUL',5X,'AUG',5X,'SEP',5X,'OCT',5X,'NOV',5X,'DEC',5X,'ANN'/)  
 2023 FORMAT (2X,'Lake St. Clair Monthly Mean Outflows',//
     &,'YEAR',3X,'JAN',5X,'FEB',5X,'MAR',5X,'APR',5X,'MAY',5X,'JUN',  
     &5X,'JUL',5X,'AUG',5X,'SEP',5X,'OCT',5X,'NOV',5X,'DEC',5X,'ANN'/)  
 2031 FORMAT (2X,'Lake Erie Monthly Mean Lake Levels',//
     &,'YEAR',3X,'JAN',5X,'FEB',5X,'MAR',5X,'APR',5X,'MAY',5X,'JUN',  
     &5X,'JUL',5X,'AUG',5X,'SEP',5X,'OCT',5X,'NOV',5X,'DEC',5X,'ANN'/)  
 2032 FORMAT (2X,'Lake Erie Beginning-of-Month Lake Levels',//
     &,'YEAR',3X,'JAN',5X,'FEB',5X,'MAR',5X,'APR',5X,'MAY',5X,'JUN',  
     &5X,'JUL',5X,'AUG',5X,'SEP',5X,'OCT',5X,'NOV',5X,'DEC',5X,'ANN'/)  
 2033 FORMAT (2X,'Lake Erie Monthly Mean Outflows - Niagara',//
     &,'YEAR',3X,'JAN',5X,'FEB',5X,'MAR',5X,'APR',5X,'MAY',5X,'JUN',  
     &5X,'JUL',5X,'AUG',5X,'SEP',5X,'OCT',5X,'NOV',5X,'DEC',5X,'ANN'/)  
 
      RETURN
      END

C-----------------------------------------------------------------------
      SUBROUTINE MONOUTLINE
     &   (unum, fm, yrstr, monv, annv, yr, m1, m2)

      IMPLICIT   NONE
      INtEGER    unum, fm, yr, m1, m2, i
      REAL       monv(12,100), annv(100), mv(12), av, t
      CHARACTER*4   yrstr

      DO i = 1, 12
         t = monv(i,yr)
         IF (i .Lt. m1) t = -10000000.
         IF (i .Gt. m2) t = -10000000.
         mv(I) = t
      ENDDO
      av = annv(yr)

      CALL WRITEMONLINE(unum, fm, yrstr, mv, av)

      END

C-----------------------------------------------------------------------
C  this subroutine prepares a single statistics line (mean, stddev, min, max)
C  for output.  It fills the monthly value array with the appropriate
C  values so that the output routine will fill the correct spots with
C  blanks (if appropriate).
C
C  there are 12 monthly values and an annual value on each line
C  the arrangement of the values is based on the start/end dates.  It could
C  be one of the 2 following choices:
C    1) _M_A    i.e. continuous set of monthly values
C    2) M_MA    i.e. blanks in middle of set  (e.g. 1-3 and 8-12)
C   where m=1 or more monthly values
C         _=0 or more blank values    (in case 2 it is 1 or more)
C         A=the annual value
C
C  the output routine will test the values, and any of them which are less
C  than -999999 will be assumed "missing" and blanks will be written in that
C  spot.
C
      SUBROUTINE STATOUTLINE
     &   (unum, fm, s, monv, annv)

      IMPLICIT   NONE
      INTEGER    unum, fm, iyr1, mon1, iyrn, monn
      INTEGER    i, totmonths
      REAL       monv(12), annv
      CHARACTER*4   s

      COMMON /PERIODS/  iyr1,mon1,iyrn,monn

      totmonths = (iyrn-iyr1)*12 + monn-mon1

C
C  If there are 12 or more months, then we have a statistical value for
C  each month of the year.
C
      IF (totmonths .GE. 12) THEN
         CALL WRITEMONLINE(unum, fm, s, monv, annv)
         RETURN
      ENDIF

C
C  Since there are less than 12 months, then one or more of the monthly
C  values is missing.  Now we need to determine which of them that is.
C  the first thing we will do is determine which layout we have.
C     i.e. is it _MA, M_A, _M_A (all handled the same)
C       or is it M_MA           (handled differently)
C  Once that is determined (a simple matter, really), then we assign the
C  "invalid" months to a value of -10000000. so that they will be written
C  as blanks by the output procedure.
C
      IF (mon1 .GT. monn) THEN     ! the layout is M_MA
         DO i = monn+1, mon1-1
           monv(i) = -10000000.
         ENDDO
      ELSE
         DO i = 1, mon1-1
           monv(i) = -10000000.
         ENDDO
         DO i = monn+1, 12
           monv(i) = -10000000.
         ENDDO
      ENDIF

      CALL WRITEMONLINE(unum, fm, s, monv, annv)

      END


C-----------------------------------------------------------------------
C  this subroutine writes out a single line into the file with 'unum' as
C  the unit number.  the format is:
C     ssssnnnnnnnnNNNNNNNNnnnnnnnnNNNNNNNN..... (A4, 13F8.x)
C  where ssss=some string such as '1987' or 'mean', 'stdv', etc
C        nnnnnnnn or NNNNNNNN = a real value in either F8.0 or F8.2 format
C  the format is based on the 'fm' variable.
C  there are 12 monthly values and an annual value on each line
C
C  We will test the values, and any of them which are less than
C  -999999 will be assumed "missing" and blanks will be written in that
C  spot.
C
      SUBROUTINE WRITEMONLINE
     &   (unum, fm, s, monv, annv)

      IMPLICIT   NONE
      INTEGER    unum, fm
      INTEGER    i,j,k
      REAL       monv(12), annv
      CHARACTER*4   s
      CHARACTER*6   fmtstr
      CHARACTER*108 line
      
      WRITE(line(1:4), 1001) s

      IF (fm .EQ. 0) fmtstr = '(F8.0)'
      IF (fm .EQ. 2) fmtstr = '(F8.2)'

      j = 5
      DO i = 1, 12
         k = J + 7
         IF (monv(i) .LT. -999999.0) THEN
            WRITE(line(j:k), '(A8)') '        '
         ELSE
            WRITE(line(j:k), fmtstr) monv(i)
         ENDIF
         j = k+1
      ENDDO

      k = j + 7
      WRITE(line(j:k), fmtstr) annv

      WRITE(unum, 1002) line      
      
 1001 FORMAT(A4)
 1002 FORMAT(A)

      END


C-----------------------------------------------------------------------
      SUBROUTINE MONSTATSLINE
     &   (unum, mons, v1, v2, v3, v4, v5, v6, v7, v8, v9)
      IMPLICIT  NONE
      INTEGER   unum, i
      REAL      v1, v2, v3, v4, v5, v6, v7, v8, v9
      CHARACTER*9  mons
      CHARACTER*78 line

      DO i=1,78
         line(i:i) = ' '
      ENDDO

      line(1:9) = mons

      IF (v1 .Gt. -999999) WRITE(line(10:16), 1001) v1
      IF (v2 .Gt. -999999) WRITE(line(18:24), 1001) v2
      IF (v3 .Gt. -999999) WRITE(line(26:32), 1002) v3
      IF (v4 .Gt. -999999) WRITE(line(33:39), 1001) v4
      IF (v5 .Gt. -999999) WRITE(line(41:47), 1001) v5
      IF (v6 .Gt. -999999) WRITE(line(49:55), 1002) v6
      IF (v7 .Gt. -999999) WRITE(line(56:62), 1001) v7
      IF (v8 .Gt. -999999) WRITE(line(64:70), 1001) v8
      IF (v9 .Gt. -999999) WRITE(line(72:78), 1002) v9

      WRITE(unum, 1010) line

 1001 FORMAT(F7.2)
 1002 FORMAT(F7.0)
 1010 FORMAT(A78)

      END

C-----------------------------------------------------------------------
      Subroutine INITDATA
C     -------------------      
C
C  Initdata reads in all input files and defines all
C  user-determined variables for a run.  Input data files
C  are all opened.

      IMPLICIT NONE
      CHARACTER title(3)*80,ext*4,run*3,type*4
      CHARACTER fnames(34)*40,units*7
      INTEGER iyr1,iyrn,mon1,monn,i,nruns,inc,nfile,mon,nyr,qm,numfile
      REAL mhk,mhym,mhsar,misar,husar,scsar,ersar
      REAL mha,sca,era,mhb,scb,erb,mhwt,scwt,erwt
      REAL sck,scym,erk,erym,dt
      REAL mhz,scz,erz,niz

C  Common block declarations.
                  
      COMMON /BOP/      mhz,scz,erz,niz
      COMMON /FILENAME/ fnames,nfile,nruns,run
      COMMON /MISC/     title,ext,type,dt,mon,nyr,inc,qm,units
      COMMON /PERIODS/  iyr1,mon1,iyrn,monn
      COMMON /SA/       mhsar,misar,husar,scsar,ersar
      COMMON /STAGE/    mhk,sck,erk,mhym,scym,erym,mha,sca,era,
     &                  mhb,scb,erb,mhwt,scwt,erwt

C  Read header and number of runs.
      
      READ(7,'(A80)') (title(I),I=1,2)
      READ(7,'(I3)') nruns

C  Read run file extension.
    
    1 READ(7,'(A4)') ext

C  Read title card.
      
      READ(7,'(A80)') title(3)

C  Read type of input data..."NBS" or "COMP"
      
      READ(7,'(/,A4)') type

C  Read initial parameters: beginning year and month, ending year and
C  month.
      
      READ(7,'(/,I4,1X,I2,1X,I4,1X,I2)') iyr1,mon1,iyrn,monn

C  Read lake surface areas.

      READ(7,'(/,5(1X,E9.4E2))') mhsar,misar,husar,scsar,ersar

C  Read units ('metric' or 'English').

      READ(7,'(/,A7)') units

C  Read discharge equation constants, outlet bottom elevations, depth and
C  fall exponents, and weighting factors for each connecting channel. 
      
      READ(7,'(/,F9.4,1X,F6.2,3(1X,F4.2))')mhk,mhym,mha,mhb,mhwt
      READ(7,'(F9.4,1X,F6.2,3(1X,F4.2))')sck,scym,sca,scb,scwt
      READ(7,'(F9.4,1X,F6.2,3(1X,F4.2))')erk,erym,era,erb,erwt

C  Read number of increments per quarter-monthly time step.
      
      READ(7,'(/,I2)') inc

C  Read initial beginning-of-month lake levels.
      
      READ (7,'(/,4F7.2,/)') mhz,scz,erz,niz

C  Determine the number of input data files.
  
      IF (type .NE. 'COMP') THEN
         numfile = 20
      ELSE
         numfile = 29
      ENDIF

C  One less input data file if Chippawa Grass Island Pool is not used.

      IF (erb .EQ. 0) THEN
         numfile = numfile-1
      ENDIF

C  Read input data file names.

      READ (7,'(A40)') (fnames(I),I=1,numfile)
      RETURN
      END

C----------------------------------------------------------------------
      Subroutine QREAD (nyear, fname, qdat)
C     ------------------------------------
C
C  Reads quarter-monthly data from the file given by fname into the array dat

      INTEGER month,nyear,year,qmon
      CHARACTER*40 fname
      REAL qdat(4,12)

      OPEN (UNIT=10, FILE=fname, STATUS='OLD', ERR=999)

C  Read data from file.

      READ (10,'(//)')
      year = -9999

      DO WHILE (year .NE. nyear)
        DO qmon = 1,4
           READ(10,'(I4,12F8.0)',ERR=999)year,(qdat(qmon,month),
     &           month=1,12)
        ENDDO
      ENDDO

      CLOSE(10)
      RETURN

 999  CALL HALT (' Error in reading file '//fname)
      CLOSE (10)

      STOP
      END

C-----------------------------------------------------------------------
      Subroutine NODAYS(iy,im,id,ndm,ndy,ierr)
C     -----------------------------------------
C
C  Nodays computes the number of days in the month
C  and the year.
C
C   Uses DATESEQUENCE routine for date calculations      
C
      IMPLICIT  NONE
      INTEGER   iy, im, id, ndm, ndy, ierr
      INTEGER   m2, y2, seq1, seq2

      ierr = 0

      IF (im .GT. 12) THEN
         WRITE (8,700)
         ierr = 1
         RETURN
      ENDIF

      IF (im .LT. 1) THEN
         WRITE (8,705)
         ierr = 1
         RETURN
      ENDIF

      IF (id .LT. 0) THEN
         WRITE (8,710)
         ierr = 1
         RETURN
      ENDIF

      IF (iy .LT. 0) THEN
         WRITE (8,715)
         ierr = 1
         RETURN
      ENDIF

c
c  How many days in the specified month?
c
      CALL DATESEQUENCE(1, im, iy, seq1)
      m2 = im + 1
      y2 = iy
      IF (m2 .GT. 12) THEN
         m2 = m2 - 12
         y2 = y2 + 1
      ENDIF
      CALL DATESEQUENCE(1, m2, y2, seq2)

      ndm = seq2 - seq1

c
c  How many days in the specified year?
c
      CALL DATESEQUENCE(1, 1, iy, seq1)
      y2 = iy + 1
      CALL DATESEQUENCE(1, 1, y2, seq2)
      ndy = seq2 - seq1

      RETURN

  700 FORMAT(1X,'MONTH GREATER THAN 12')
  705 FORMAT(1X,'MONTH IS TOO SMALL')
  710 FORMAT (1X,'DAY IS NEGATIVE')
  715 FORMAT (1X,'YEAR IS NEGATIVE')
  720 FORMAT (1X,'DAY IS TOO LARGE')

      END


C-------------------------------------------------------------------------
      SUBROUTINE DATESEQUENCE
     &  (day, month, year, sequencenumber)
C
C     This routine computes the number of the date, with days numbered
C     sequentially from 1 (corresponding to the base date).  the base
C     date must be 1 March on any LEAP YEAR or century boundary. Dates
C     must be between the base date and the base date plus 2,147,483,646.
C     This routine is set up to use only 4-byte integer variables and 
C     to avoid intermediate use of reals in the evaluation of expressions.
C     It is a modification of the previous version which used 2-byte integer
C     variables. For reference on the 2-byte version of this subroutine
C     see Dr. Dobb's Journal 80:66-70.  Adapted by Tim Hunter.
C
C     Note that this routine correctly handles the year 2000, as well as all
C     other leap year issues.
C
      IMPLICIT NONE
      INTEGER   baseyear
      PARAMETER (baseyear = 0)
      INTEGER   numyears, totc, leapc, skipc, nonleapc
      INTEGER   day, month, year, sequencenumber
      IF (year .LE. baseyear) STOP 'year .LE. baseyear'
      numyears = year - baseyear - 1 + (month + 9) / 12
      totc = (baseyear + numyears)/100 - (baseyear/100)
      skipc = 3 - MOD(baseyear/100, 4)
      leapc = (totc + 3 - skipc) / 4
      nonleapc = totc - leapc
      sequencenumber = 365 * numyears + (numyears/4) - nonleapc
     &                + (153 * MOD(month + 9, 12) + 2) / 5 + day 
      RETURN
      END

C-----------------------------------------------------------------------
      Subroutine CLRSCR ()
C     --------------------
C
C  CLRSCR clears the terminal screen.

      IMPLICIT LOGICAL (a - Z)
      INTEGER loop

      DO 10 loop = 1, 25
         PRINT *
   10 CONTINUE

      RETURN
      END

C-----------------------------------------------------------------------
      Subroutine HALT(message)
C     ------------------------
C
C  HALT clears the screen, prints the given error message, and stops
C  execution.


      CHARACTER message*80

      CALL CLRSCR

      WRITE(*,'(1X,/)')
      WRITE(*,'(1x,A80)') message

      STOP
      END



