!=============================================================================
!----------------------------------------------------------------------------
!     LARGE BASIN RUNOFF MODEL (6 June 2000)
!
!     Contact Thomas E. Croley II or Timothy S. Hunter
!             (734) 741-2238         (734) 741-2344
!             (734) 741-2055 (FAX)   (734) 741-2055 (FAX)
!             croley@glerl.noaa.gov  hunter@glerl.noaa.gov
!
!             Great Lakes Environmental Research Laboratory
!             National Oceanic and Atmospheric Administration
!             US Department of Commerce
!             2205 Commonwealth Blvd.
!             Ann Arbor, Michigan  48105-2945
!
!     This program is the general-purpose version of the Large Basin
!     Runoff Model (LBRM) as developed by Dr. Thomas E. Croley II.  The
!     original Fortran 77 model code, as specifically adapted for use in
!     GLERL's Advanced Hydrologic Prediction System (AHPS), has been
!     adapted by Tim Hunter (Dec 1999) to be more generic.  It was also
!     translated to Fortran 95 structure in order to take advantage of
!     MODULEs and eliminate COMMON blocks.  This version of the LBRM is
!     distributed along with calibration programs and documentation
!     [Croley, T.E., II, 2000.  Large basin runoff model.  Chapter in
!     "Mathematical Models of Watershed Hydrology" (D. Frevert, S. Meyer,
!     and V.P. Singh, eds.), Water Resources Publications, Littleton,
!     Colorado (in prep.).]  All earlier versions of the LBRM were
!     specific to GLERL uses, so this program was developed.  This
!     version is designed so that it will work for any arbitrary
!     watershed (not tied to the Great Lakes subbasins).
!
!     SPECIAL NOTE:
!     This program code is designed to be used in two ways: the first is as a
!     Windows DLL to be called by a Windows program; the second is as a stand-
!     alone DOS program.  Because of the absence of a standard set of condition-
!     al compiler directives, the user will have to compile and link in two dif-
!     ferent ways, to build either a DLL or a DOS program.  Furthermore, the
!     user will have to compile and link differently for DOS depending on use of
!     the program in a simulation setting or a calibration setting.
!
!     The DOS commands to compile this code and link it for the three cases are as
!     follows (using the Lahey fortran 95 compiler as an example):
!
!     First Case (for use as Windows DLL):
!
!     lf95 ErrorProcessWIN.f95 -c -nco -win -dll
!     lf95 GLERLutil.f95       -c -nco -win -dll
!     lf95 LBRM.f95            -c -nco -win -dll
!     lf95 LBRM_BD.f95         -c -nco -win -dll -ml bd
!     lf95 MyMsgAPI.f95        -c -nco -win -dll -ml winapi
!     lf95 MyMsgAPI.obj LBRM_BD.obj LBRM.obj GLERLutil.obj ErrorProcessWIN -win -dll -out CALIBRUN.DLL
!
!     Second Case (for use in SIMULATION as a DOS program):
!
!     lf95 ErrorProcessDOS.f95 -c
!     lf95 GLERLutil.f95       -c
!     lf95 LBRM.f95            -c
!     lf95 DOSmodel.f95        -c
!     lf95 DOSmodel.obj LBRM.obj GLERLutil.obj ErrorProcessDOS -out DOSmodel.EXE
!
!     Third Case (for use in CALIBRATION as a DOS program):
!
!     lf95 ErrorProcessDOS.f95 -c
!     lf95 GLERLutil.f95       -c
!     lf95 LBRM.f95            -c
!     lf95 DOScalib.f95        -c
!     lf95 DOScalib.obj LBRM.obj GLERLutil.obj ErrorProcessDOS -out DOScalib.EXE
!
!
!     FILE FORMAT DESCRIPTIONS
!     ------------------------
!
!     LBRM.CFG
!        This is the configuration (or control) file which contains the
!        file names.  That's all it has in it.  The format is similar to
!        the older Windows ".INI" files.  That is, all lines will be ignored
!        except for those that start with the key string we are looking for
!        followed by an "=" sign.  The rest of the line will be considered
!        a file name.  An example file might look like the following:
!
!        |** This is the configuration file for the LBRM **
!        |This line will be ignored.
!        |$$$ So will this one.
!        |   spaces at the beginning of a line are insignificant.
!        | The ordering of the lines is also insignificant.
!        |The next 4 lines are the only ones we need to have in here.
!        |PARMFILE=PARAM.TXT
!        |BCFILE = BNDRCOND.TXT
!        |  METFILE =   METDATA.TXT
!        |OUTFILE = OUTPUT.TXT
!        |DATAFILE =FLOWDATA.TXT
!        | WORKDIR =  C:\LBRM\
!
!
!        Note that the entry "WORKDIR=" was included for use by the Windows-based
!        applications.  It is not used in this Fortran code and will be ignored.
!
!     PARAMETER FILE  (could be named anything)
!        This file contains what would be considered the control information
!        for the LBRM in simulation or calibration.  That includes dates,
!        watershed area, and parameters.
!
!        The format is as follows:
!
!        1st record: Identifier (FORMAT: A)
!                    |        |
!                    |________|___Any string the user wants to use which will
!                                 help identify this model run. This will be
!                                 echoed into the output file simply for the
!                                 user's help later.  It can be left blank with
!                                 no consequences if the user so desires.
!
!        2nd record: Area (FORMAT: E13.6E2)
!                    |  |
!                    |__|____watershed area in square meters.
!
!        3rd record: iD iM iY iD iM iY (FORMAT: 2I2, I4, 2I2, I4)
!                    || || || || || ||
!                    || || || || || ||___starting year for calibration run
!                    || || || || ||______starting month for calibration run
!                    || || || ||_________starting day for calibration run
!                    || || ||
!                    || || ||   NOTE: calibration date unused in simulation
!                    || || ||
!                    || || ||___starting year in data file for model run
!                    || ||______starting month in data file for model run
!                    ||_________starting day in data file for model run
!
!        4th record: iD iM iY (FORMAT: 2I2, I4)
!                    || || ||
!                    || || ||___ending year in data file for model run
!                    || ||______ending month in data file for model run
!                    ||_________ending day in data file for model run
!
!        5th through 13th records: model parameter and calibration info
!
!                    Parameter I J Minimum Maximum (FORMAT: E13.6E2)
!                    |       | | | |     | |     |
!                    |       | | | |     | |_____|_allowable parameter maximum (calibration only)
!                    |       | | | |_____|_allowable parameter minimum (calibration only)
!                    |       | | |_maximum number significant digits (calibration only)
!                    |       | |_starting significant digit (calibration only)
!                    |_______|_model parameter as follows:
!
!             Tbase   base temperature for computation of heat available
!                     for evapotranspiration
!             AlbedS  melt factor (proportionality constant) for computation
!                     of degree-day snowmelt
!             AlpPer  linear reservoir coefficient for percolation from
!                     the upper soil zone moisture storage
!             AlpUEv  partial linear reservoir coefficient for evapo-
!                     transpiration from the upper soil zone
!             AlpInt  linear reservoir coefficient for interflow from the
!                     lower soil zone moisture storage
!             AlpDPr  linear reservoir coefficient for deep percolation
!                     from the lower soil zone moisture storage
!             AlpLEv  partial linear reservoir coefficient for evapo-
!                     transpiration from the lower soil zone
!             AlpGw   linear reservoir coefficient for groundwater flow
!                     from the groundwater zone moisture storage
!             AlpSf   linear reservoir coefficient for basin outflow from
!                     the surface zone moisture storage
!
!        14th through 29th records: (FORMAT: E13.6E2)
!
!             Cons    proportionality constant for computation of heat
!                     available for evapotranspiration; CANNOT BE ZERO DURING
!                     PRODUCTION RUNS; IF ZERO DURING CALIBRATION, "Cons" IS
!                     COMPUTED FROM THE FOLLOWING 14 VALUES:
!             Rjan    January average mid-monthly daily cloudless solar
! 		      insolation at the surface (langleys)
!             Rfeb    February insolation
!             Rmar    March insolation
!             Rapr    April insolation
!             Rmay    May insolation
!             Rjun    June insolation
!             Rjul    July insolation
!             Raug    August insolation
!             Rsep    September insolation
!             Roct    October insolation
!             Rnov    November insolation
!             Rdec    December insolation
!             B1      empirical linear coefficient [ordinate-intercept]
!		      relating ratio of cloudy insolation/cloudless insolation
!                     to "cloud cover"
!             B2      empirical linear coefficient [slope] relating ratio of
!                     cloudy insolation/cloudless insolation to "cloud cover"
!             USZC    upper soil zone capacity (cm)
!
!       30th record: iParm Corr RMSE (FORMAT: I2, 2E13.6E2)
!                    |   | |  | |  |
!                    |   | |  | |__|__root mean square error (calibration only)
!                    |   | |__|_______correlation (calibration only)
!                    |___|____________parameter in start calibration with (calibration only)
!
!     BOUNDARY CONDITIONS FILE (could be named anything)
!        This file contains boundary conditions that the model outputs must adhere
!        to.  During the model run, values of the various storages are (re)set to
!        the boundary conditions on the indicated dates if the boundary conditions
!        are non-negative.  A date, corresponding to the starting date in the
!        parameter file, must be present in this file and would, by definition,
!        contain the initial conditions (which CANNOT be negative).  There may be
!        any number of boundary conditions present in this file, but THEY MUST BE
!        IN CHRONOLOGICAL ORDER!
!
!        1st record: skipped (may be anything)
!
!        2nd record: skipped (may be anything)
!
!        3-? records: iD iM iY Uszm Lszm Gzm Ss Snw (FORMAT: 2I2, I4, 5E13.6E2)
!                     || || || |  | |  | | | || | |
!                     || || || |  | |  | | | || |_|___snowpack moisture (cm.)
!                     || || || |  | |  | | | ||_______surface storage (cm.)
!                     || || || |  | |  | |_|__________groundwater storage (cm.)
!                     || || || |  | |__|______________lower soil zone storage (cm.)
!                     || || || |__|___________________upper soil zone storage (cm.)
!                     || || ||________________________date's year
!                     || ||___________________________date's month
!                     ||______________________________date's day
!
!
!     METEOROLOGICAL DATA FILE  (could be named anything)
!        This file contains the daily air temperatures (maximum and minimum)
!        and precipitation for each day of the model run period specified
!        in the parameter file.  It may contain data for dates outside that
!        period as well.  If so, only the data for the specified period will
!        actually be used.
!
!        Note that NO MISSING DATA IS ALLOWED!!!!!!  If the LBRM encounters
!        missing data during a run it will immediately stop with an error message.
!
!        The format is as follows:
!
!        1st record: Identifier (FORMAT: A)
!                    |        |
!                    |________|___Any string the user wants to use which will
!                                 help identify this data set.  It will not be
!                                 used by the LBRM for anything.
!
!        2nd record: "DDMMYYYY TMIN TMAX PREC" (FORMAT: A)
!                     |                     |
!                     |_____________________|__header info string - not used
!
!        Records 3 through the end of the file have daily data.  These records
!        must be in chronological order with no missing days.
!
!        dd mm yyyy Tmin Tmax Precip (FORMAT: 2I2, I4, 3I5)
!        || || |  | |  | |  | |    |
!        || || |  | |  | |  | |____|__depth of daily precipitation in
!        || || |  | |  | |  |         hundredths of millimeters.
!        || || |  | |  | |__|_________maximum daily air temperature in 
!        || || |  | |  |              hundredths of degrees Celsius.
!        || || |  | |__|______________minimum daily air temperature in
!        || || |  |                   hundredths of degrees Celsius.
!        || || |__|___________________four-digit year of data's date
!        || ||________________________two-digit month of data's date
!        ||___________________________two-digit day of data's date
!
!
!     OUTPUT DATA FILE  (could be named anything)
!        This file contains the daily output from the LBRM for each day
!        of the model run period specified in the parameter file.
!
!        The format is as follows:
!
!        1st record: Identifier (FORMAT: A)
!                    |        |
!                    |________|___This is the string that was read from record
!                                 1 of the PARAMETER file.
!
!        2nd record: "Output from the LBRM.  Runoff and storages are in mm over the basin."
!
!        3rd record: "DDMMYYYY  RUNOFF   USZM   LSZM GROUND SURFAC   SNOW"
!
!        Records 4 through the end of the file have the daily output (end-of-day).
!
!           dd mm yyyy Run USZ LSZ GZM SRF SNW (FORMAT: 2I2, I4, 6E11.4E2)
!           || || |  | | | | | | | | | | | | |
!           || || |  | | | | | | | | | | | |_|__snowpack moisture storage (cm)
!           || || |  | | | | | | | | | |_|______surface moisture storage (cm)
!           || || |  | | | | | | | |_|__________groundwater moisture storage (cm)
!           || || |  | | | | | |_|______________lower soil zone moisture storage (cm)
!           || || |  | | | |_|__________________upper soil zone moisture storage (cm)
!           || || |  | |_|______________________runoff volume for this day (cm)
!           || || |__|__________________________four-digit year of data's date
!           || ||_______________________________two-digit month of data's date
!           ||__________________________________two-digit day of data's date
!
!-----------------------------------------------------------------
MODULE LBRM_Global

   USE GlerlUtil
   
!
!  If this is being used in GLERL AHPS package then we will use file paths as
!  specified in the AHPS configuration file.  If this program is being used in
!  a "stand-alone" application then we set this logical to FALSE and will not
!  use the file name info.  The settings of these variables will be left up to
!  the main calling programs in either case.
!
   LOGICAL            :: UseAhpsCfg
   CHARACTER (LEN=99) :: AhpsCfgFile

END MODULE LBRM_Global

!===========================================================================
!===========================================================================
!=============================================================================

MODULE LBRM_Subs

   USE LBRM_Global
!
!     In Fortran 77 source code, these were in an unnamed COMMON block
!
!     COMMON Uszm, Lszm, Gzm, Ss, Evap, Hplse, UszmAvg, LszmAvg, USZC,
!    & AlpPer, AlpUEv, AlpInt, AlpDPr, AlpLEv, AlpGw, AlpGEv, AlpSf,
!    & AlpSEv, CID, Epsilon, Dpsilon, Gpsilon, Hpsilon, vRun, vInt,
!    & vPer, vGw, vUEv, vLEv
!
!     AlpDPr = linear reservoir constant for deep percolation, inv. days
!     AlpGEv = partial constant of groundwater evaporation, inv. cub. m.
!     AlpGw  = linear reservoir constant for grooundwater flow, inv. days
!     AlpInt = linear reservoir constant for interflow, inv. days
!     AlpLEv = partail constant of lower zone evaporation, inv. cub. m.
!     AlpPer = linear reservoir constant for percolation, inv. days
!     AlpSEv = partial constant of surface evaporation, inv. cub. m.
!     AlpSf  = linear reservoir constant for surface flow, inv. days
!     AlpUEv = partial constant of upper zone evap., inv. cub. m.
!     CID    = time in one group of days (week, month, etc.), days
!     Evap   = total evapotranspiration volume, cub. m.
!     Gzm    = groundwater zone moisture, cub. m.
!     Hplse  = total energy out (Evap. + Pot. Evap.) water equ., cub. m.
!     Lszm   = lower soil zone moisture, cub. m.
!     LszmAvg= average lower soil zone moisture storage over the time period.
!     Ss     = surface water storage, cub. m.
!     USZC   = upper soil zone moisture capacity, cub. m.
!     Uszm   = upper soil zone moisture, cub. m.
!     UszmAvg= average upper soil zone moisture storage over the time period.
!     vGw    = groundwater zone outflow volume, cub. m.
!     vInt   = interflow volume, cub. m.
!     vLEv   = lower zone evapotranspiration volume, cub. m.
!     vPer   = percolation volume, cub. m.
!     vRun   = surface runoff volume, cub. m.
!     vUEv   = upper zone evapotranspiration volume, cub. m.
!
!     Epsilon, Dpsilon, Gpsilon, and Hpsilon : convergence criteria
!
      REAL  :: AlpDPr, AlpGEv, AlpGw, AlpInt, AlpLEv, AlpPer
      REAL  :: AlpSEv, AlpSf, AlpUEv
      REAL  :: CID, Evap, Gzm, Hplse, Lszm, LszmAvg
      REAL  :: Ss, USZC, Uszm, UszmAvg
      REAL  :: vGw, vInt, vLEv, vPer, vRun, vUEv
      REAL  :: Epsilon, Dpsilon, Gpsilon, Hpsilon

CONTAINS
!-----------------------------------------------------------------
      SUBROUTINE OUTFLOW (NS)
!
!     This routine computes the mass balance for the GLERL Large Basin Runoff
!     Model, a tank cascade model of basin storages and runoff.
!
      IMPLICIT NONE
      REAL, INTENT(INOUT) :: NS
      REAL  :: EvPRp, NSR, vDPr, vGEv, vInf
      REAL  :: A, B, C, D, E, F, G, H, L, M, N, O, P, R, T
!
!      COMMON Uszm, Lszm, Gzm, Ss, Evap, Hplse, UszmAvg, LszmAvg, USZC,
!     & AlpPer, AlpUEv, AlpInt, AlpDPr, AlpLEv, AlpGw, AlpGEv, AlpSf,
!     & AlpSEv, CID, Epsilon, Dpsilon, Gpsilon, Hpsilon, vRun, vInt,
!     & vPer, vGw, vUEv, vLEv
!
!     EvPRp = potential evapotranspiration rate, cub. m./ day
!     NS    = on input  : net supply volume to the watershed surface, cub. m.
!             on output : basin outflow volume, cub. m. on output
!     NSR   = net supply rate, cub. m. / day
!     vDPr  = deep percolation volume, cub. m.
!     vGEv  = groundwater zone evapotranspiration volume, cub. m.
!     vInf  = infiltration volume, cub. m.
!
!     A, B, C, D, E, F, G, H, L, M,   \ /  dummy variables for storage
!     N, O, P, R, and T               / \  of intermediate results
!
      NSR = NS / CID
      EvPRp = Hplse / 2. / CID
!
!     Use iterative subroutine to determine EvPRp, if necessary.
!
      IF (EvPRp .GT. Dpsilon) CALL EvpRate(NSR, EvPRp)
      B = NSR/USZC + AlpPer + AlpUEv * EvPRp
      C = AlpPer * NSR / B
      A = AlpPer * Uszm - C
      T = EPN(-B * CID) * A / AlpPer + NSR / B
      UszmAvg = (A / B * (1. - EPN(-B * CID)) / CID + C) / AlpPer
      R = NS + Uszm - T
      Uszm = T
      vInf = NS - R * NSR / USZC / B
      vPer = R * AlpPer / B
      vRun = NS - vInf
      vUEv = R - vRun - vPer
      D = AlpInt + AlpDPr + Alplev * EvPRp
      IF (ABS((D - B) / D) .LE. Epsilon .OR. ABS(D - B) .LE. Hpsilon) THEN
        T = (Lszm + A * CID - C / D) * EPN(-D * CID) + C / D
        F = A
        G = C / D
        E = Lszm - G
        LszmAvg = (E / D * (1. - EPN(-D * CID)) + F /D ** 2 * (1. -       &
     &            (D * CID + 1.) * EPN(-D * CID))) / CID + G
        R = vPer + Lszm - T
        IF (AlpDPr .GE. Dpsilon) THEN
          E = AlpDPr * E
          F = AlpDPr * F
          G = AlpDPr * G
        ENDIF
        Lszm = T
        vInt = R * AlpInt / D
        vDPr = R * AlpDPr / D
        vLEv = R - vInt - vDPr
        H = AlpGw + AlpGEv * EvPRp
        IF (ABS((H - D) / H) .LE. Epsilon .OR. ABS(H - D) .LE. Hpsilon) THEN
          T = (Gzm + E * CID + F/2. * CID**2 - G/H) * EPN(-H * CID) + G/H
          IF (AlpDPr .LT. Dpsilon) T = Gzm * EPN(-H * CID)
          R = vDPr + Gzm - T
          IF (AlpDPr .LT. Dpsilon) THEN
            L = AlpInt * E + A * NSR / USZC / AlpPer
            M = AlpInt * F
            N = 0.
            O = AlpInt * G + C * NSR / USZC / AlpPer
          ELSE
            L = AlpGw * (Gzm - G / H) + AlpInt / AlpDPr * E +             &
     &          NSR / USZC / AlpPer * A
            M = AlpGw * E + AlpInt / AlpDPr * F
            N = AlpGw * F / 2.
            O = C * NSR / USZC / AlpPer + G * (AlpInt / AlpDPr + AlpGw / H)
          ENDIF
          Gzm = T
          vGw = R * AlpGw / H
          vGEv = R - vGw
          P = AlpSf + AlpSEv * EvPRp
          IF (ABS((P-H) / P) .LE. Epsilon .OR. ABS(P-H) .LE. Hpsilon) THEN
            T = (Ss + L * CID + M / 2. * CID ** 2 + N / 3. * CID ** 3     &
     &          - O / P) * EPN(-P * CID) + O / P
          ELSE
            T = (Ss - L / (P-H) + M / (P-H) ** 2 - 2. * N / (P-H)**3 -    &
     &          O/P) * EPN(-P * CID) + (L / (P-H) + M / (P-H)**2 *        &
     &          ((P-H) * CID - 1.) + N / (P-H) * CID**2 - 2. * N /        &
     &          (P-H)**3 * ((P - H) * CID - 1.)) * EPN(-H * CID) + O / P
          ENDIF
        ELSE
          T = (Gzm - E / (H-D) + F / (H-D)**2 - G/H) * EPN(-H * CID) +    &
     &        (E / (H-D) + F / (H-D) ** 2 * ((H-D) * CID - 1.)) *         &
     &        EPN(-D * CID) + G / H
          IF (AlpDPr .LT. Dpsilon) T = Gzm * EPN(-H * CID)
          R = vDPr + Gzm - T
          IF (AlpDPr .LT. Dpsilon) THEN
            L = 0.
            M = AlpInt * E + A * NSR / USZC / AlpPer
            N = AlpInt * F
            O = AlpInt * G + C * NSR / USZC / AlpPer
          ELSE
            L = AlpGw * (Gzm + F / (H - D) ** 2 - E / (H - D) - G / H)
            M = AlpGw * (E / (H - D) - F / (H - D) ** 2) + AlpInt         &
     &          / AlpDPr * E + (NSR / USZC / AlpPer) * A
            N = F * (AlpGw / (H - D) + AlpInt / AlpDPr)
            O = C * NSR / USZC / AlpPer + G * (AlpInt/AlpDPr + AlpGw/H)
          ENDIF
          Gzm = T
          vGw = R * AlpGw / H
          vGEv = R - vGw
          P = AlpSf + AlpSEv * EvPRp
          IF (ABS((P-H) / P) .LE. Epsilon .OR. ABS(P-H) .LE. Hpsilon) THEN
            T = (Ss + L * CID - M / (P-D) + N / (P-D) ** 2 - O / P) *     &
     &          EPN(-P * CID) + (M / (P-D) + N / (P-D)**2 *               &
     &          ((P-D) * CID - 1.)) * EPN(-D * CID) + O / P
          ELSE
            IF (ABS((P-D) / P) .LE. Epsilon .OR. ABS(P-D) .LE. Hpsilon) THEN
              T = (Ss - L / (P-H) + M * CID + N / 2. * CID**2 - O/P) *    &
     &            EPN(-P * CID) + L / (P - H) * EPN(-H * CID) + O / P
            ELSE
              T = (Ss - L / (P-H) - M / (P-D) + N / (P-D) ** 2 - O/P) *   &
     &            EPN(-P * CID) + L / (P - H) * EPN(-H * CID) + O / P +   &
     &            (M/(P-D) + N/(P-D)**2 * ((P-D)*CID - 1.)) * EPN(-D*CID)
            ENDIF
          ENDIF
        ENDIF
      ELSE
        F = A / (D - B)
        G = C / D
        E = Lszm - F - G
        T = E * EPN(-D * CID) + F * EPN(-B * CID) + G
        LszmAvg = (E / D * (1. - EPN(-D * CID)) + F / B *                 &
     &            (1. - EPN(-B * CID))) / CID + G
        R = vPer + Lszm - T
        IF (AlpDPr .GE. Dpsilon) THEN
          E = AlpDPr * E
          F = AlpDPr * F
          G = AlpDPr * G
        ENDIF
        Lszm = T
        vInt = R * AlpInt / D
        vDPr = R * AlpDPr / D
        vLEv = R - vInt - vDPr
        H = AlpGw + AlpGEv * EvPRp
        IF (ABS((H-D) / H) .LE. Epsilon .OR. ABS(H-D) .LE. Hpsilon) THEN
          T = (Gzm + E * CID - F / (H - B) - G / H) * EPN(-H * CID) +     &
     &        F / (H - B) * EPN(-B * CID) + G / H
          IF (AlpDPr .LT. Dpsilon) T = Gzm * EPN(-H * CID)
          R = vDPr + Gzm - T
          IF (AlpDPr .LT. Dpsilon) THEN
            L = AlpInt * E
            M = 0.
            N = AlpInt * F + A * NSR / USZC / AlpPer
            O = AlpInt * G + C * NSR / USZC / AlpPer
          ELSE
            L = AlpGw * (Gzm - F / (H - B) - G / H) + AlpInt / AlpDPr * E
            M = AlpGw * E
            N = A * NSR/USZC/AlpPer + F * (AlpInt/AlpDPr + AlpGw / (H-B))
            O = C * NSR/USZC/AlpPer + G * (AlpInt/AlpDPr + AlpGw / H)
          ENDIF
          Gzm = T
          vGw = R * AlpGw / H
          vGEv = R - vGw
          P = AlpSf + AlpSEv * EvPRp
          IF (ABS((P-H) / P) .LE. Epsilon .OR. ABS(P-H) .LE. Hpsilon) THEN
            T = (Ss + L * CID + M / 2. * CID**2 - N / (P-B) - O / P) *    &
     &          EPN(-P * CID) + N / (P - B) * EPN(-B * CID) + O / P
          ELSE
            IF (ABS((P-B) / P) .LE. Epsilon .OR. ABS(P-B) .LE. Hpsilon) THEN
              T = (Ss - L / (P-H) + M / (P-H)**2 + N * CID - O/P) *       &
     &            EPN(-P * CID) + (L / (P-H) + M / (P-H)**2 *             &
     &            ((P-H) * CID - 1.)) * EPN(-H * CID) + O / P
            ELSE
              T = (Ss - L / (P-H) + M / (P-H)**2 - N / (P-B) - O/P) *     &
     &            EPN(-P * CID) + (L / (P-H) + M / (P-H)**2 *             &
     &            ((P-H) * CID - 1.)) * EPN(-H * CID) + N / (P - B) *     &
     &            EPN(-B * CID) + O / P
            ENDIF
          ENDIF
        ELSE
          IF (ABS((H-B) / H) .LE. Epsilon .OR. ABS(H-B) .LE. Hpsilon) THEN
            T = (Gzm - E / (H-D) + F * CID - G / H) * EPN(-H * CID) +     &
     &          E / (H - D) * EPN(-D * CID) + G / H
            IF (AlpDPr .LT. Dpsilon) T = Gzm * EPN(-H * CID)
            R = vDPr + Gzm - T
            IF (AlpDPr .LT. Dpsilon) THEN
              L = AlpInt * F + A * NSR / USZC / AlpPer
              M = AlpInt * E
              N = 0.
              O = AlpInt * G + C * NSR / USZC / AlpPer
            ELSE
              L = AlpGw * (Gzm - E / (H-D) - G / H) +                     &
     &            AlpInt / AlpDPr * F + A * NSR / USZC / AlpPer
              M = E * (AlpInt / AlpDPr + AlpGw / (H - D))
              N = AlpGw * F
              O = C * NSR/USZC/AlpPer + G * (AlpInt/AlpDPr + AlpGw/H)
            ENDIF
            Gzm = T
            vGw = R * AlpGw / H
            vGEv = R - vGw
            P = AlpSf + AlpSEv * EvPRp
            IF (ABS((P-H) / P) .LE. Epsilon .OR. ABS(P-H) .LE. Hpsilon) THEN
              T = (Ss + L * CID - M / (P-D) + N / 2. * CID**2 - O/P) *    &
     &            EPN(-P * CID) + M / (P-D) * EPN(-D * CID) + O/P
            ELSE
              IF (ABS((P-D)/P) .LE. Epsilon .OR. ABS(P-D) .LE. Hpsilon) THEN
                T = (Ss - L / (P-H) + M * CID + N / (P-H)**2 - O/P) *     &
     &              EPN(-P * CID) + (L / (P-H) + N / (P-H)**2  *          &
     &              ((P-H) * CID - 1.)) * EPN(-H * CID) + O/P
              ELSE
                T = (Ss - L / (P-H) - M / (P-D) + N / (P-H)**2 - O/P) *   &
     &              EPN(-P * CID) + (L / (P-H) + N / (P-H)**2 *           &
     &              ((P-H) * CID - 1.)) * EPN(-H * CID) + M / (P-D) *     &
     &              EPN(-D * CID) + O / P
              ENDIF
            ENDIF
          ELSE
            T = (Gzm - E / (H-D) - F / (H-B) - G/H) * EPN(-H * CID) +     &
     &          E / (H-D) * EPN(-D * CID) + F/(H-B) * EPN(-B * CID) + G/H
            IF (AlpDPr .LT. Dpsilon) T = Gzm * EPN(-H * CID)
            R = vDPr + Gzm - T
            IF (AlpDPr .LT. Dpsilon) THEN
              L = 0.
              M = AlpInt * E
              N = AlpInt * F + A * NSR / USZC / AlpPer
              O = AlpInt * G + C * NSR / USZC / AlpPer
            ELSE
              L = AlpGw * (Gzm - E / (H - D) - F / (H - B) - G / H)
              M = E * (AlpInt / AlpDPr + AlpGw / (H - D))
              N = A * NSR/USZC/AlpPer + F * (AlpInt/AlpDPr + AlpGw/(H-B))
              O = C * NSR/USZC/AlpPer + G * (AlpInt/AlpDPr + AlpGw/H)
            ENDIF
            Gzm = T
            vGw = R * AlpGw / H
            vGEv = R - vGw
            P = AlpSf + AlpSEv * EvPRp
            IF (ABS((P-H) / P) .LE. Epsilon .OR. ABS(P-H) .LE. Hpsilon) THEN
              T = (Ss + L * CID - M / (P-D) - N / (P-B) - O/P) *          &
     &            EPN(-P * CID) + M / (P-D) * EPN(-D * CID) + N/(P-B) *   &
     &            EPN(-B * CID) + O / P
            ELSE
              IF (ABS((P-D)/P) .LE. Epsilon .OR. ABS(P-D) .LE. Hpsilon) THEN
                T = (Ss - L / (P-H) + M * CID - N / (P-B) - O/P) *        &
     &              EPN(-P * CID) + L / (P-H) * EPN(-H * CID) + N /       &
     &              (P-B) * EPN(-B * CID) + O / P
              ELSE
                IF (ABS((P-B)/P) .LE. Epsilon .OR. ABS(P-B) .LE. Hpsilon) THEN
                  T = (Ss - L / (P-H) - M / (P-D) + N*CID - O/P) *        &
     &                EPN(-P * CID) + L / (P-H) * EPN(-H * CID) +         &
     &                M / (P-D) * EPN(-D * CID) + O / P
                ELSE
                  T = (Ss - L / (P-H) - M / (P-D) - N / (P-B) - O/P) *    &
     &                EPN(-P * CID) + L / (P-H) * EPN(-H * CID) + M /     &
     &                (P-D) * EPN(-D*CID) + N/(P-B) * EPN(-B*CID) + O/P
                ENDIF
              ENDIF
            ENDIF
          ENDIF
        ENDIF
      ENDIF
      R = vRun + vInt + vGw + Ss - T
      Ns = R * AlpSf / P
      Ss = T
      Evap = vUEv + vLEv + vGEv + R - Ns
      RETURN
      END SUBROUTINE OUTFLOW
!-------------------------------------------------------------------------
      SUBROUTINE EvpRate(NSR, Evp)
      IMPLICIT NONE
      REAL, INTENT(INOUT)  ::  NSR, Evp
      REAL (KIND=REAL8) :: A, B, C, D, E, F, G
      REAL (KIND=REAL8) :: uzAvg, lzAvg, EvPRp, EvPRpo
!
!      COMMON Uszm, Lszm, Gzm, Ss, Evap, Hplse, UszmAvg, LszmAvg, USZC,
!     & AlpPer, AlpUEv, AlpInt, AlpDPr, AlpLEv, AlpGw, AlpGEv, AlpSf,
!     & AlpSEv, CID, Epsilon, Dpsilon, Gpsilon, Hpsilon, vRun, vInt,
!     & vPer, vGw, vUEv, vLEv
!
!     Determine EvPRp via iteration.
!
      EvPRp = Evp
      EvPRpo = 0.
      DO WHILE (ABS((EvPRp - EvPRpo) / EvPRp) .GT. Gpsilon)
        B = NSR / USZC + AlpPer + AlpUEv * EvPRp
        C = NSR / B
        A = Uszm - C
        uzAvg = A / B * (1. - EPX(-B * CID)) / CID + C
        D = AlpInt + AlpDPr + AlpLEv * EvPRp
        IF (ABS((D-B)/D) .LE. Epsilon .OR. ABS(D-B) .LE. Hpsilon) THEN
          F = AlpPer * A
          G = AlpPer * C / D
          E = Lszm - G
          lzAvg = (E / D * (1. - EPX(-D * CID)) + F / D**2 *              &
     &            (1. - (D * CID + 1.) * EPX(-D * CID))) / CID + G
        ELSE
          F = AlpPer * A / (D - B)
          G = AlpPer * C / D
          E = Lszm - F - G
          lzAvg = (E / D * (1. - EPX(-D * CID)) + F / B *                 &
     &            (1. - EPX(-B * CID))) / CID + G
        ENDIF
        EvPRpo = EvPRp
        EvPRp = Hplse / CID / (1. + AlpUEv * uzAvg + AlpLEv * lzAvg)
      END DO
      Evp = EvPRp
      RETURN
      END SUBROUTINE EvpRate
!--------------------------------------------------------------------------
      REAL FUNCTION EPN(X)  ! designed for use with negative arguments
      IMPLICIT NONE
      REAL, INTENT(IN) :: X
      IF (X .LT. -87.33) THEN
        EPN = 0.
      ELSE
!        IF (X .GT. 88.72) THEN
!           EPN = HUGE(X)
!        ELSE
           EPN = EXP(X)
!        END IF
      ENDIF
      RETURN
      END FUNCTION EPN
!--------------------------------------------------------------------------
      REAL (KIND=REAL8) FUNCTION EPX(X)  ! designed for use with negative arguments
      IMPLICIT NONE
      REAL (KIND=REAL8), INTENT(IN) :: X
      IF (X .LT. -708.39) THEN
         EPX = 0.
      ELSE
 !        IF (X .GT. 709.77) THEN
 !           EPX = HUGE(X)
 !        ELSE
            EPX = EXP(X)
 !        END IF
      END IF
      RETURN
      END FUNCTION EPX
!--------------------------------------------------------------------------
END MODULE LBRM_Subs

!===========================================================================
!===========================================================================
!=============================================================================
MODULE LBRM_IO

   USE LBRM_Global

      INTEGER :: iNoBC, iNoDataPoints
      INTEGER :: iFDMR, iLDMR, iFDCR, iParm, iNoChng, NumBC, iNoDays
      INTEGER, DIMENSION(9) :: Ixrc
      REAL :: Area, SvRMSE, CorrS
      REAL, DIMENSION(25) :: Parameter
      CHARACTER (LEN=200) :: CnfgFile, ParmFile, MetFile, OutFile, BCfile, DataFile
      CHARACTER (LEN=99) :: RunID
      CHARACTER (LEN=180), DIMENSION(29) :: Cmnt

CONTAINS
!--------------------------------------------------------------------------
      SUBROUTINE FindArraySizes(RunMetSize, BoundSize)
      IMPLICIT NONE
      INTEGER, INTENT(INOUT) :: RunMetSize, BoundSize
      INTEGER :: IOS, iD, iM, iY, Seq, SSeq, ESeq, U1, U2
      CHARACTER (LEN=200) :: Line

      ErrorLevel = 0
!
!     Compute number of entries to be put into the runoff (output)
!     data file.  This is the same as the number of entries to be
!     read from the meteorological data (input) file.
!
      RunMetSize = 0
      U1 = GetFreeUnitNumber(); IF (ErrorLevel .NE. 0) GOTO 899
      OPEN(UNIT=U1, FILE=TRIM(ParmFile), STATUS='OLD', ERR=811)
      CALL FileWasOpened(U1); IF (ErrorLevel .NE. 0) GOTO 899
      READ(U1, *, ERR=812)
      READ(U1, *, ERR=812)
      READ(U1, 1005, ERR=812) iD, iM, iY
      CALL DateSequence(iD, iM, iY, SSeq); IF (ErrorLevel .NE. 0) GOTO 899
      READ(U1, 1005, ERR=812) iD, iM, iY
      CALL DateSequence(iD, iM, iY, ESeq); IF (ErrorLevel .NE. 0) GOTO 899
      CLOSE(U1)
      CALL FileWasClosed(U1); IF (ErrorLevel .NE. 0) GOTO 899
      RunMetSize = ESeq - SSeq + 1
!
!     Find number of entries in the boundary conditions file.
!     Only count those between the desired simulation start and end dates
!     that were just read from the parameter file (SSeq & ESeq).
!
      BoundSize = 0
      U2 = GetFreeUnitNumber(); IF (ErrorLevel .NE. 0) GOTO 899
      OPEN(UNIT=U2, FILE=TRIM(BCFile), STATUS='OLD', ERR=821)
      CALL FileWasOpened(U2); IF (ErrorLevel .NE. 0) GOTO 899
      READ(U2, *, ERR=822)
      READ(U2, *, ERR=822)
      READ(U2, 1001, IOSTAT=IOS) Line
      DO WHILE (IOS .EQ. 0)
         READ(Line, 1005, ERR=899) iD, iM, iY
         CALL DateSequence(iD, iM, iY, Seq); IF (ErrorLevel .NE. 0) GOTO 899
         IF (Seq .GE. SSeq .AND. Seq .LE. ESeq) BoundSize = BoundSize + 1 
         READ(U2, 1001, IOSTAT=IOS) Line
      END DO
      CLOSE(U2)
      CALL FileWasClosed(U2); IF (ErrorLevel .NE. 0) GOTO 899
      RETURN

  811 ErrorMessage = 'LBRM: Error opening file '//TRIM(ParmFile)
      CALL PassMsg; GOTO 898
  812 ErrorMessage = 'LBRM: Error reading file '//TRIM(ParmFile)
      CALL PassMsg; GOTO 898
  821 ErrorMessage = 'LBRM: Error opening file '//TRIM(BCFile)
      CALL PassMsg; GOTO 898
  822 ErrorMessage = 'LBRM: Error reading file '//TRIM(BCFile)
      CALL PassMsg; GOTO 898

  898 ErrorLevel = 1  ! Return with an error indication.
  899 RETURN

 1001 FORMAT(A)
 1005 FORMAT(2I2, I4)
      END SUBROUTINE FindArraySizes
!--------------------------------------------------------------------------
      SUBROUTINE ReadLBRMCFGfile
      IMPLICIT NONE
      INTEGER :: IOS, U1
      CHARACTER (LEN=10)  :: Item
      CHARACTER (LEN=180) :: Value
      CHARACTER (LEN=200) :: Line

      ErrorLevel = 0

!
!     Set file path variable for config file if needed.  All this requires
!     is that if the global flag variable has been set to FALSE by the 
!     calling procedure (main program) then force the string variable to be
!     completely empty.  In this case the assumption will be that the LBRM.CFG
!     file is in the current directory and that the main prog already put
!     us in the right directory via CHDIR commands or whatever.
!
      IF (.NOT. UseAhpsCfg) AhpsCfgFile = ''

!
!     Read information from the LBRM.CFG file
!
      ParmFile = '---'
      MetFile  = '---'
      OutFile  = '---'
      BCfile   = '---'
      DataFile = '---'
      CnfgFile = TRIM(AhpsCfgFile) // 'LBRM.CFG'

      U1 = GetFreeUnitNumber(); IF (ErrorLevel .NE. 0) GOTO 899
      OPEN(UNIT=U1, FILE=TRIM(CnfgFile), STATUS='OLD', ERR=811)
      CALL FileWasOpened(U1); IF (ErrorLevel .NE. 0) GOTO 899
      READ(U1, 1101, IOSTAT=IOS) Line
      DO WHILE (IOS .EQ. 0)
         CALL Caps(Line); IF (ErrorLevel .NE. 0) GOTO 899
         CALL ParseCfgLine(Line, Item, Value); IF (Errorlevel .NE. 0) GOTO 899
         IF (Item(1:8) .EQ. 'PARMFILE') ParmFile = TRIM(Value)
         IF (Item(1:7) .EQ. 'METFILE')  MetFile = TRIM(Value)
         IF (Item(1:7) .EQ. 'OUTFILE')  OutFile = TRIM(Value)
         IF (Item(1:6) .EQ. 'BCFILE')   BCfile = TRIM(Value)
         IF (Item(1:8) .EQ. 'DATAFILE') DataFile = TRIM(Value)
         READ(U1, 1101, IOSTAT=IOS) Line
      END DO
      IF (IOS .GT. 0) GOTO 812     ! error other than eof
      CLOSE(U1)
      CALL FileWasClosed(U1); IF (ErrorLevel .NE. 0) GOTO 899

      IF (ParmFile(1:3) .EQ. '---') THEN
         ErrorMessage = 'No PARMFILE entry in LBRM.CFG!'
         CALL PassMsg; GOTO 898
      ENDIF
      IF (MetFile(1:3) .EQ. '---') THEN
         ErrorMessage = 'No METFILE entry in LBRM.CFG!'
         CALL PassMsg; GOTO 898
      ENDIF
      IF (OutFile(1:3) .EQ. '---') THEN
         ErrorMessage = 'No OUTFILE entry in LBRM.CFG!'
         CALL PassMsg; GOTO 898
      ENDIF
      IF (BCfile(1:3) .EQ. '---') THEN
         ErrorMessage = 'No BCFILE entry in LBRM.CFG!'
         CALL PassMsg; GOTO 898
      ENDIF
      RETURN

  811 ErrorMessage = 'LBRM: Error opening file '//TRIM(CnfgFile)
      CALL PassMsg; GOTO 898
  812 ErrorMessage = 'LBRM: Error reading file '//TRIM(CnfgFile)
      CALL PassMsg; GOTO 898
  898 ErrorLevel = 1         ! Return with an error indication.
  899 RETURN

 1101 FORMAT(A)
      END SUBROUTINE ReadLBRMCFGfile
!--------------------------------------------------------------------------
      SUBROUTINE ReadParameterFile
      IMPLICIT NONE
      INTEGER :: iD, iM, iY, I1, I2, I3, U1

      ErrorLevel = 0
!
!     Read values from the parameter file
!
      U1 = GetFreeUnitNumber(); IF (ErrorLevel .NE. 0) GOTO 899
      OPEN(UNIT=U1, FILE=TRIM(ParmFile), STATUS='OLD', ERR=811)
      CALL FileWasOpened(U1); IF (ErrorLevel .NE. 0) GOTO 899
      READ(U1, 2101, ERR=812) RunID
      READ(U1, 2102, ERR=812) Area, Cmnt(1)
!
!     Start and end dates
!
      READ(U1, 2103, ERR=812) iD, iM, iY, I1, I2, I3, Cmnt(2)
      CALL DateSequence(iD, iM, iY, iFDMR); IF (ErrorLevel .NE. 0) GOTO 899
      CALL DateSequence(I1, I2, I3, iFDCR); IF (ErrorLevel .NE. 0) GOTO 899
      READ(U1, 2104, ERR=812) iD, iM, iY, Cmnt(3)
      CALL DateSequence(iD, iM, iY, iLDMR); IF (ErrorLevel .NE. 0) GOTO 899
      iNoDays = iLDMR - iFDMR + 1
      IF (iNoDays .GT. iNoDataPoints) THEN
         ErrorMessage = 'LBRM: Simulation dates exceed available storage'
         CALL PassMsg; GOTO 898
      ENDIF
      IF (iNoDays .LT. 1) THEN
         ErrorMessage = 'LBRM: Simulation dates are invalid'
         CALL PassMsg; GOTO 898
      ENDIF
!
!     Parameter values
!
      READ(U1, 2105, ERR=812) Parameter( 1), Ixrc(1), iNoChng, Cmnt( 4)  ! TBase
      READ(U1, 2105, ERR=812) Parameter( 2), Ixrc(2), iD, Cmnt( 5)  ! AlbedS
      READ(U1, 2105, ERR=812) Parameter( 3), Ixrc(3), iD, Cmnt( 6)  ! AlpPer
      READ(U1, 2105, ERR=812) Parameter( 4), Ixrc(4), iD, Cmnt( 7)  ! AlpUEv
      READ(U1, 2105, ERR=812) Parameter( 5), Ixrc(5), iD, Cmnt( 8)  ! AlpInt
      READ(U1, 2105, ERR=812) Parameter( 6), Ixrc(6), iD, Cmnt( 9)  ! AlpDpr
      READ(U1, 2105, ERR=812) Parameter( 7), Ixrc(7), iD, Cmnt(10)  ! AlpLev
      READ(U1, 2105, ERR=812) Parameter( 8), Ixrc(8), iD, Cmnt(11)  ! AlpGw
      READ(U1, 2105, ERR=812) Parameter( 9), Ixrc(9), iD, Cmnt(12)  ! AlpSf
      READ(U1, 2102, ERR=812) Parameter(10), Cmnt(13)  ! Cons
      READ(U1, 2102, ERR=812) Parameter(11), Cmnt(14)  ! R( 2)
      READ(U1, 2102, ERR=812) Parameter(12), Cmnt(15)  ! R( 3)
      READ(U1, 2102, ERR=812) Parameter(13), Cmnt(16)  ! R( 4)
      READ(U1, 2102, ERR=812) Parameter(14), Cmnt(17)  ! R( 5)
      READ(U1, 2102, ERR=812) Parameter(15), Cmnt(18)  ! R( 6)
      READ(U1, 2102, ERR=812) Parameter(16), Cmnt(19)  ! R( 7)
      READ(U1, 2102, ERR=812) Parameter(17), Cmnt(20)  ! R( 8)
      READ(U1, 2102, ERR=812) Parameter(18), Cmnt(21)  ! R( 9)
      READ(U1, 2102, ERR=812) Parameter(19), Cmnt(22)  ! R(10)
      READ(U1, 2102, ERR=812) Parameter(20), Cmnt(23)  ! R(11)
      READ(U1, 2102, ERR=812) Parameter(21), Cmnt(24)  ! R(12)
      READ(U1, 2102, ERR=812) Parameter(22), Cmnt(25)  ! R(13)
      READ(U1, 2102, ERR=812) Parameter(23), Cmnt(26)  ! B1          ! 0.355
      READ(U1, 2102, ERR=812) Parameter(24), Cmnt(27)  ! B2          ! 0.68
      READ(U1, 2102, ERR=812) Parameter(25), Cmnt(28)  ! USZC
      READ(U1, 2106, ERR=812) iParm, SvRMSE, CorrS, Cmnt(29)
      CLOSE(U1)
      CALL FileWasClosed(U1); IF (ErrorLevel .NE. 0) GOTO 899
      RETURN

  811 ErrorMessage = 'LBRM: Error opening file '//TRIM(ParmFile)
      CALL PassMsg; GOTO 898
  812 ErrorMessage = 'LBRM: Error reading file '//TRIM(ParmFile)
      CALL PassMsg; GOTO 898

  898 ErrorLevel = 1  ! Return with an error indication.
  899 RETURN

 2101 FORMAT(A)
 2102 FORMAT(E13.6E2, A)
 2103 FORMAT(2I2, I4, 2I2, I4, A)
 2104 FORMAT(2I2, I4, A)
 2105 FORMAT(E13.6E2, 2I2, A)
 2106 FORMAT(I2, 2E13.6E2, A)
      END SUBROUTINE ReadParameterFile
!--------------------------------------------------------------------------
      SUBROUTINE WriteParameterFile
      IMPLICIT NONE
      INTEGER :: IOS, iD, iM, iY, I1, I2, I3, U1
!
!     Write values to the parameter file
!
      U1 = GetFreeUnitNumber(); IF (ErrorLevel .NE. 0) GOTO 899
      OPEN(UNIT=U1, FILE=TRIM(ParmFile), STATUS='OLD', ERR=10)
      CLOSE(UNIT=U1, STATUS='DELETE')
  10  OPEN(UNIT=U1, FILE=TRIM(ParmFile), STATUS='NEW', ERR=811)
      CALL FileWasOpened(U1); IF (ErrorLevel .NE. 0) GOTO 899

      WRITE(U1, 2101, ERR=813) TRIM(RunID)
      WRITE(U1, 2102, ERR=813) Area, TRIM(Cmnt(1))
      CALL SequenceDate(iD, iM, iY, iFDMR); IF (ErrorLevel .NE. 0) GOTO 899
      CALL SequenceDate(I1, I2, I3, iFDCR); IF (ErrorLevel .NE. 0) GOTO 899
      WRITE(U1, 2103, ERR=813) iD, iM, iY, I1, I2, I3, TRIM(Cmnt(2))
      CALL SequenceDate(iD, iM, iY, iLDMR); IF (ErrorLevel .NE. 0) GOTO 899
      WRITE(U1, 2104, ERR=813) iD, iM, iY, TRIM(Cmnt(3))
      WRITE(U1, 2105, ERR=813) Parameter( 1), Ixrc(1), iNoChng, TRIM(Cmnt( 4))
      WRITE(U1, 2105, ERR=813) Parameter( 2), Ixrc(2), iNoChng, TRIM(Cmnt( 5))
      WRITE(U1, 2105, ERR=813) Parameter( 3), Ixrc(3), iNoChng, TRIM(Cmnt( 6))
      WRITE(U1, 2105, ERR=813) Parameter( 4), Ixrc(4), iNoChng, TRIM(Cmnt( 7))
      WRITE(U1, 2105, ERR=813) Parameter( 5), Ixrc(5), iNoChng, TRIM(Cmnt( 8))
      WRITE(U1, 2105, ERR=813) Parameter( 6), Ixrc(6), iNoChng, TRIM(Cmnt( 9))
      WRITE(U1, 2105, ERR=813) Parameter( 7), Ixrc(7), iNoChng, TRIM(Cmnt(10))
      WRITE(U1, 2105, ERR=813) Parameter( 8), Ixrc(8), iNoChng, TRIM(Cmnt(11))
      WRITE(U1, 2105, ERR=813) Parameter( 9), Ixrc(9), iNoChng, TRIM(Cmnt(12))
      WRITE(U1, 2102, ERR=813) Parameter(10), TRIM(Cmnt(13))
      WRITE(U1, 2102, ERR=813) Parameter(11), TRIM(Cmnt(14))
      WRITE(U1, 2102, ERR=813) Parameter(12), TRIM(Cmnt(15))
      WRITE(U1, 2102, ERR=813) Parameter(13), TRIM(Cmnt(16))
      WRITE(U1, 2102, ERR=813) Parameter(14), TRIM(Cmnt(17))
      WRITE(U1, 2102, ERR=813) Parameter(15), TRIM(Cmnt(18))
      WRITE(U1, 2102, ERR=813) Parameter(16), TRIM(Cmnt(19))
      WRITE(U1, 2102, ERR=813) Parameter(17), TRIM(Cmnt(20))
      WRITE(U1, 2102, ERR=813) Parameter(18), TRIM(Cmnt(21))
      WRITE(U1, 2102, ERR=813) Parameter(19), TRIM(Cmnt(22))
      WRITE(U1, 2102, ERR=813) Parameter(20), TRIM(Cmnt(23))
      WRITE(U1, 2102, ERR=813) Parameter(21), TRIM(Cmnt(24))
      WRITE(U1, 2102, ERR=813) Parameter(22), TRIM(Cmnt(25))
      WRITE(U1, 2102, ERR=813) Parameter(23), TRIM(Cmnt(26))
      WRITE(U1, 2102, ERR=813) Parameter(24), TRIM(Cmnt(27))
      WRITE(U1, 2102, ERR=813) Parameter(25), TRIM(Cmnt(28))
      WRITE(U1, 2106, ERR=813) iParm, SvRMSE, CorrS, Cmnt(29)
      CLOSE(U1)
      CALL FileWasClosed(U1); IF (ErrorLevel .NE. 0) GOTO 899
      RETURN

  811 ErrorMessage = 'LBRM: Error opening file '//TRIM(ParmFile)
      CALL PassMsg; GOTO 898
  813 ErrorMessage = 'LBRM: Error writing file '//TRIM(ParmFile)
      CALL PassMsg; GOTO 898

  898 ErrorLevel = 1  ! Return with an error indication.
  899 RETURN

 2101 FORMAT(A)
 2102 FORMAT(E13.6E2, A)
 2103 FORMAT(2I2, I4, 2I2, I4, A)
 2104 FORMAT(2I2, I4, A)
 2105 FORMAT(E13.6E2, 2I2, A)
 2106 FORMAT(I2, 2E13.6E2, A)
      END SUBROUTINE WriteParameterFile
!--------------------------------------------------------------------------
      SUBROUTINE ReadBoundaryConditionsFile(iBCD, Upper, Lower, Ground, Surface, Snow)
      IMPLICIT NONE
      INTEGER, INTENT(INOUT), DIMENSION(:) :: iBCD
      REAL,    INTENT(INOUT), DIMENSION(:) :: Upper, Lower, Ground, Surface, Snow
      INTEGER :: IOS, iD, iM, iY, Seq, U1
      REAL :: Uszm, Lszm, Gzm, Ss, Snw
      CHARACTER (LEN=200) :: Line

      ErrorLevel = 0
!
!     Read values from the boundary conditions file
!
      U1 = GetFreeUnitNumber(); IF (ErrorLevel .NE. 0) GOTO 899
      OPEN(UNIT=U1, FILE=TRIM(BCfile), STATUS='OLD', ERR=811)
      CALL FileWasOpened(U1); IF (ErrorLevel .NE. 0) GOTO 899
      READ(U1, 2101, ERR=812)
      READ(U1, 2102, ERR=812)
!
!     Save any boundary conditions available on or after start date
!     and on or before end date.
!
      NumBC = 0
      READ(U1, 2101, IOSTAT=IOS) Line
      DO WHILE (IOS .EQ. 0)
         READ(Line, 2115, ERR=899) iD, iM, iY, Uszm, Lszm, Gzm, Ss, Snw
         CALL DateSequence(iD, iM, iY, Seq); IF (ErrorLevel .NE. 0) GOTO 899
         IF (Seq .GE. iFDMR .AND. Seq .LE. iLDMR) THEN
            NumBC = NumBC + 1  ! number of boundary conditions read in and saved
            IF (NumBC .GT. iNoBC) THEN
               ErrorMessage = 'Number of boundary conditions exceeds allowable'
               CALL PassMsg; GOTO 898
            ENDIF
            iBCD(NumBC) = Seq
!
!           Convert from centimeters over watershed to cubic meters (missing
!           data denoted negative values still interpreted as negative values).
!
            Upper(NumBC) = Uszm * Area / 100.
            Lower(NumBC) = Lszm * Area / 100.
            Ground(NumBC) = Gzm * Area / 100.
            Surface(NumBC) = Ss * Area / 100.
            Snow(NumBC) = Snw * Area / 100.
         ENDIF
         READ(U1, 2101, IOSTAT=IOS) Line
      END DO
      CLOSE(U1)
      CALL FileWasClosed(U1); IF (ErrorLevel .NE. 0) GOTO 899
      IF (IOS .GT. 0) GOTO 812
      RETURN

  811 ErrorMessage = 'LBRM: Error opening file '//TRIM(BCFile)
      CALL PassMsg; GOTO 898
  812 ErrorMessage = 'LBRM: Error reading file '//TRIM(BCFile)
      CALL PassMsg; GOTO 898

  898 ErrorLevel = 1  ! Return with an error indication.
  899 RETURN

 2101 FORMAT(A)
 2102 FORMAT(E13.6E2)
 2115 FORMAT(2I2, I4, 5E13.6E2)
      END SUBROUTINE ReadBoundaryConditionsFile
!--------------------------------------------------------------------------
      SUBROUTINE ReadLBRMDataFile(Tmx, Tmn, Prc)
      IMPLICIT NONE
      INTEGER (KIND=INTEGER2), INTENT(INOUT), DIMENSION(:) :: Tmx, Tmn, Prc
      INTEGER :: I, IOS, iD, iM, iY, Seq, I1, I2, I3, U1

      ErrorLevel = 0
!
!     Read all of the meteorological data.  We only need it one day at 
!     a time, and the original program read it as needed in order to
!     conserve precious memory resources under the DOS 640K limit.  We 
!     are no longer constrained by that limit, though, and the additional
!     memory that is required to store the met data for the entire period
!     (2 bytes/value * 3 data types = 6 bytes/day) is not a problem.
!     100 years of data only requires 214K.
!
!     By reading the data in "one fell swoop" we ensure that hard disk 
!     "head thrashing" doesn't occur.  It also means that we can do all 
!     of the error checking on the input data now rather than finding a
!     problem after we've already processed 90% of the period.
!
      U1 = GetFreeUnitNumber(); IF (ErrorLevel .NE. 0) GOTO 899
      OPEN(UNIT=U1, FILE=TRIM(MetFile), STATUS='OLD', ERR=811)
      CALL FileWasOpened(U1); IF (ErrorLevel .NE. 0) GOTO 899
      READ(U1, *, ERR=812)
      READ(U1, *, ERR=812)
      READ(U1, 3103, ERR=812) iD, iM, iY, I1, I2, I3
      CALL DateSequence(iD, iM, iY, Seq); IF (ErrorLevel .NE. 0) GOTO 899
      IF (Seq .GT. iFDMR) THEN
         ErrorMessage = TRIM(MetFile)//' does not sufficiently contain date range.'
         CALL PassMsg;  GOTO 898
      END IF
      IOS = 0
      I = Seq - iFDMR
      DO WHILE (IOS .EQ. 0 .AND. I .LT. iNoDays)
         CALL DateSequence(iD, iM, iY, Seq); IF (ErrorLevel .NE. 0) GOTO 899
         IF (Seq - iFDMR .NE. I) THEN
            ErrorMessage = TRIM(MetFile)//' is not chronological or is discontinuous.'
            CALL PassMsg; GOTO 898
         END IF
         I = I + 1
         IF (I .GE. 1) THEN
            Tmn(I) = I1
            Tmx(I) = I2
            Prc(I) = I3
         END IF
         READ(U1, 3103, IOSTAT=IOS) iD, iM, iY, I1, I2, I3
         IF (IOS .GT. 0) GOTO 812
      END DO
      CLOSE(U1)
      CALL FileWasClosed(U1); IF (ErrorLevel .NE. 0) GOTO 899
      IF (I .NE. iNoDays) THEN
         ErrorMessage = TRIM(MetFile)//' does not sufficiently contain date range.'
         CALL PassMsg;  GOTO 898
      END IF
      DO I = 1, iNoDays
        IF ((Tmn(I).EQ.-9999) .OR. (Tmx(I).EQ.-9999) .OR. (Prc(I).EQ.-9999)) THEN
          ErrorMessage = 'LBRM: Missing data in '//TRIM(MetFile)//' is not allowed!'
          CALL PassMsg;  GOTO 898
        END IF
      END DO
      RETURN

  811 ErrorMessage = 'LBRM: Error opening file '//TRIM(MetFile)
      CALL PassMsg; GOTO 898
  812 ErrorMessage = 'LBRM: Error reading file '//TRIM(MetFile)
      CALL PassMsg; GOTO 898

  898 ErrorLevel = 1  ! Return with an error indication.
  899 RETURN

 3103 FORMAT(2I2, I4, 3I5)
      END SUBROUTINE ReadLBRMDataFile
!--------------------------------------------------------------------------
      SUBROUTINE InputLBRMdata(BC, DS, iBCD, Upper, Lower, Ground,        &
     &                         Surface, Snow, Tmx, TMn, Prc)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: BC, DS
      INTEGER (KIND=INTEGER2), DIMENSION(DS), INTENT(INOUT) :: Tmx, Tmn, Prc
      INTEGER, DIMENSION(BC), INTENT(INOUT) :: iBCD
      REAL,    DIMENSION(BC), INTENT(INOUT) :: Upper,Lower,Ground,Surface,Snow
!
!     Read information from the LBRM.CFG file
!
      CALL ReadLBRMCFGfile
      IF (ErrorLevel .NE. 0) RETURN
!
!     Read values from the parameter file
!
      CALL ReadParameterFile
      IF (ErrorLevel .NE. 0) RETURN
!
!     Read values from the boundary conditions file.
!
      CALL ReadBoundaryConditionsFile(iBCD, Upper, Lower, Ground, Surface, Snow)
      IF (ErrorLevel .NE. 0) RETURN
!
!     Read all of the meteorological data.
!
      CALL ReadLBRMDataFile(Tmx, Tmn, Prc)
      RETURN
      END SUBROUTINE InputLBRMData
!--------------------------------------------------------------------------
      SUBROUTINE WriteLBRMOutputFile(OutRunoff, OutUszm, OutLszm,       &
     &                               OutGzm, OutSs, OutSnw)
      IMPLICIT NONE
      REAL, INTENT(IN), DIMENSION(:) :: OutRunoff, OutUszm, OutLszm
      REAL, INTENT(IN), DIMENSION(:) :: OutGzm, OutSs, OutSnw
      INTEGER :: I, IOS, iD, iM, iY, Seq, U1
      REAL :: Runoff, Uszm, Lszm, Gzm, Ss, Snw

      ErrorLevel = 0
!
!     Open the file for output
!
      U1 = GetFreeUnitNumber(); IF (ErrorLevel .NE. 0) GOTO 899
      OPEN(UNIT=U1, FILE=TRIM(OutFile), STATUS='OLD', ERR=10)
      CLOSE(U1, STATUS='DELETE')
 10   OPEN(UNIT=U1, FILE=TRIM(OutFile), STATUS='NEW', ERR=811)
      CALL FileWasOpened(U1); IF (ErrorLevel .NE. 0) GOTO 899
      WRITE(U1, 4101, ERR=813) RunID
      WRITE(U1, 4102, ERR=813)
      WRITE(U1, 4103, ERR=813)
      DO I = 1, iNoDays
         Seq = I + iFDMR - 1
         CALL SequenceDate(iD, iM, iY, Seq); IF (ErrorLevel .NE. 0) GOTO 899
!
!        The runoff and storages are converted from cubic meters to
!        centimeters over the watershed area.
!
         Runoff = OutRunoff(I) / Area * 100.
         Uszm   = OutUszm(I)   / Area * 100.
         Lszm   = OutLszm(I)   / Area * 100.
         Gzm    = OutGzm(I)    / Area * 100.
         Ss     = OutSs(I)     / Area * 100.
         Snw    = OutSnw(I)    / Area * 100.
         IF (Runoff .LT. 1.2e-38) Runoff = 0.
         IF (Uszm   .LT. 1.2e-38) Uszm   = 0.
         IF (Lszm   .LT. 1.2e-38) Lszm   = 0.
         IF (Gzm    .LT. 1.2e-38) Gzm    = 0.
         IF (Ss     .LT. 1.2e-38) Ss     = 0.
         IF (Snw    .LT. 1.2e-38) Snw    = 0.
         WRITE(U1, 4104, ERR=813) iD, iM, iY, Runoff, Uszm, Lszm, Gzm, Ss, Snw
      END DO
      CLOSE (U1)
      CALL FileWasClosed(U1);  IF (ErrorLevel .NE. 0) GOTO 899
      RETURN

  811 ErrorMessage = 'LBRM: Error opening file '//TRIM(OutFile)
      CALL PassMsg; GOTO 898
  813 ErrorMessage = 'LBRM: Error writing file '//TRIM(OutFile)
      CALL PassMsg; GOTO 898

  898 ErrorLevel = 1  ! Return with an error indication.
  899 RETURN

 4101 FORMAT(A)
 4102 FORMAT('Output from the LBRM.  Runoff and storages are in cm over the basin.')
 4103 FORMAT('DDMMYYYY     RUNOFF       USZM       LSZM     GROUND    SURFACE       SNOW')
 4104 FORMAT(2I2, I4, 6E11.4E2)
      END SUBROUTINE WriteLBRMOutputFile
!--------------------------------------------------------------------------
END MODULE LBRM_IO

!===========================================================================
!===========================================================================
!===========================================================================
MODULE ArraysForLBRM

     USE LBRM_IO

     INTEGER :: BC, DS
!
!     These will be dimensioned to the number of boundary conditions (BC)
!
      INTEGER, DIMENSION(:), ALLOCATABLE :: iBCD
      REAL, DIMENSION(:), ALLOCATABLE :: Upper, Lower, Ground, Surface, Snow
!
!     These will be dimensioned to the number of days in the simulation (DS)
!
      INTEGER (KIND=INTEGER2), DIMENSION(:), ALLOCATABLE :: Tmx, Tmn, Prc
      REAL, DIMENSION(:), ALLOCATABLE :: OutRunoff, OutUszm, OutLszm
      REAL, DIMENSION(:), ALLOCATABLE :: OutGzm, OutSs, OutSnw

CONTAINS
!--------------------------------------------------------------------------
      SUBROUTINE AllocateFromDOS
      INTEGER  :: IOS

      CALL ReadLBRMCFGfile           ! Read information from the LBRM.CFG file
      IF (ErrorLevel .NE. 0) RETURN
      CALL FindArraySizes(DS, BC)    ! Find sizes based on start/end dates
      IF (ErrorLevel .NE. 0) RETURN
      iNoDataPoints = DS
      iNoBC = BC
!
!     Actually allocate the RAM
!
      ALLOCATE(iBCD(BC), Upper(BC), Lower(BC), Ground(BC),                &
     &         Surface(BC), Snow(BC), Tmx(DS), TMn(DS), Prc(DS),          &
     &         OutRunoff(DS), OutUszm(DS), OutLszm(DS), OutGzm(DS),       &
     &         OutSs(DS), OutSnw(DS), STAT=IOS)
      IF (IOS .NE. 0) THEN
         ErrorMessage = 'LBRM: Error allocating RAM for arrays.'
         CALL PassMsg
         ErrorLevel = 1
      ENDIF
      RETURN
      END SUBROUTINE AllocateFromDOS
!-------------------------------------------------
      SUBROUTINE DeallocateFromDOS
      INTEGER  :: IOS

      DEALLOCATE(OutRunoff, OutUszm, OutLszm, OutGzm, OutSs, OutSnw,      &
     &      Tmx, TMn, Prc, iBCD, Upper, Lower, Ground, Surface, Snow,     &
     &      STAT=IOS)
      END SUBROUTINE DeallocateFromDOS
!-------------------------------------------------
END MODULE ArraysForLBRM

!===========================================================================
!===========================================================================
!===========================================================================
MODULE LBRM_Main

   USE ArraysForLBRM
   USE LBRM_Subs
!   USE LBRM_IO

CONTAINS
!--------------------------------------------------------------------------
!
!       Program for interactive investigation of model
!
!
!     Area   = watershed area, sq. m.
!     iD     = calendar day of the year
!     iFDMR  = first date for model run (ddmmyyyy)
!     iLDMR  = last date for model run (ddmmyyyy)
!     iM     = calendar month of the year
!     iNDPGOD= number of days per group of days, days
!     iNoDays= number of days to be considered in data set, days
!     iY     = calendar year
!     Precip = daily precipitation volume (liquid equivalent), cub. m.
!     Snw    = snowpack volume (liquid equivalent), cub. m.
!     Ta     = average daily air temperature, deg. C.
!     Tmax   = maximum daily air temperature, deg. C.
!     Tmin   = minimum daily air temperature, deg. C.
!
!     COMMON Uszm, Lszm, Gzm, Ss, Evap, Hplse, UszmAvg, LszmAvg, USZC,
!    & AlpPer, AlpUEv, AlpInt, AlpDPr, AlpLEv, AlpGw, AlpGEv, AlpSf,
!    & AlpSEv, CID, Epsilon, Dpsilon, Gpsilon, Hpsilon, vRun, vInt,
!    & vPer, vGw, vUEv, vLEv
!
      SUBROUTINE LBRM_Engine(cCons, BC, DS,                               &
     &                       iBCD, Upper, Lower, Ground, Surface, Snow,   &
     &                       Tmx, TMn, Prc,                               &
     &                       OutRunoff, OutUszm, OutLszm, OutGzm, OutSs,  &
     &                       OutSnw,                                      &
     &                       LTnsg, LTrun, LTint, LTgwf, LTuev, LTlev)
!
!     If "cCons" is a zero, then this program is being used in a
!     "production run;" i.e., it is being used with an already-calibrated
!     parameter file in a simulation setting.  If "cCons" is nonzero,
!     then this program is being used in a "calibration run;" i.e.,
!     it is being called with alternate values of the parameter file in an
!     investigation of the effects of different parameter values on the
!     model output.
!
!     In a production run, the value of "Cons," that is read from the parameter
!     file, is NEVER allowed to be zero!  In a calibration run, if "Cons" is
!     read in as zero, then it is recomputed internally from other values in the
!     parameter file and from data in the meteorological data file.  The value
!     of "Cons" (in both cases) is returned in the above subroutine header as
!     "cCons."
!
!     LTnsg, LTrun, LTint, LTgwf, LTuev, and LTlev are respectively long-term
!     volumes of net supply to the ground surface (snowmelt and/or precipitation),
!     surface runoff, interflow, groundwater flow, upper zone evapotranspiration,
!     and lower zone evapotranspiration.
!
      IMPLICIT NONE
      INTEGER, INTENT(INOUT) :: BC, DS
      REAL,    INTENT(INOUT) :: cCons
      INTEGER (KIND=INTEGER2), DIMENSION(DS), INTENT(INOUT) :: Tmx, Tmn, Prc
      INTEGER, DIMENSION(BC), INTENT(INOUT) :: iBCD
      REAL,    DIMENSION(BC), INTENT(INOUT) :: Upper,Lower,Ground,Surface,Snow
      REAL, DIMENSION(DS), INTENT(INOUT) :: OutRunoff, OutUszm, OutLszm
      REAL, DIMENSION(DS), INTENT(INOUT) :: OutGzm, OutSs, OutSnw
      REAL, INTENT(INOUT) :: LTnsg, LTrun, LTint, LTgwf, LTuev, LTlev
      INTEGER :: I, J, iD, iM, iY, iNDPGOD
      INTEGER :: Today, II
      REAL :: Precip, Snw, Ta, Tmax, Tmin, Tbase, AlbedS
      REAL :: Cons, Melt, Runoff
      REAL :: B1, B2, AvgRR, AvgHPlE, RR
      REAL, DIMENSION(14) :: R

      ErrorLevel = 0
!
!     Set parameters from cached inputs.
!
      TBase  = Parameter( 1)
      AlbedS = Parameter( 2) * Area / 100.
      AlpPer = Parameter( 3)
      AlpUEv = Parameter( 4)
      AlpInt = Parameter( 5)
      AlpDpr = Parameter( 6)
      AlpLev = Parameter( 7)
      AlpGw  = Parameter( 8)
      AlpSf  = Parameter( 9)
      Cons   = Parameter(10) * Area * 10000.
      R( 2)  = Parameter(11)
      R( 3)  = Parameter(12)
      R( 4)  = Parameter(13)
      R( 5)  = Parameter(14)
      R( 6)  = Parameter(15)
      R( 7)  = Parameter(16)
      R( 8)  = Parameter(17)
      R( 9)  = Parameter(18)
      R(10)  = Parameter(19)
      R(11)  = Parameter(20)
      R(12)  = Parameter(21)
      R(13)  = Parameter(22)
      B1     = Parameter(23)
      B2     = Parameter(24)
      USZC   = Parameter(25)
      R( 1) = R(13)
      R(14) = R( 2)
!
!     Set boundary (and initial) conditions from cached inputs.
!
      IF ((NumBC .LE. 0) .OR. (iBCD(1) .NE. iFDMR)) THEN
         ErrorMessage = 'Initial Conditions are not present on start date'//     &
     &                  ' in ' // TRIM(BCFile)
         CALL PassMsg; GOTO 898
      END IF
      Uszm = Upper(1)
      Lszm = Lower(1)
      Gzm  = Ground(1)
      Ss   = Surface(1)
      Snw  = Snow(1)
      IF ((Uszm .LT. 0) .OR. (Lszm .LT. 0) .OR. (Gzm .LT. 0) .OR.         &
     &    (Ss   .LT. 0) .OR. (Snw  .LT. 0)) THEN
         ErrorMessage = 'Initial Conditions are not present on start date'//     &
     &                  ' in ' // TRIM(BCFile)
         CALL PassMsg; GOTO 898
      END IF

!
!     Set initial values for some variables
!
      Epsilon = 1.e-5
      Dpsilon = 1.e-30
      Gpsilon = 1.e-3
      Hpsilon = 1.e-6
      AlpGEv = 0.
      AlpSEv = 0.
!
!     Always uses 1-Day computation interval.
!
      iNDPGOD = 1
      CID = FLOAT(iNDPGOD)      ! Computational Interval Days
!     USZC   = Area * 0.02
      USZC = USZC * Area / 100.
!
!     Check and compute, if necessary, a value for "Cons"
!
      IF (cCons .EQ. 0) THEN         ! in a "production" run
         IF (Cons .EQ. 0) THEN       ! not allowed in a production run
            ErrorMessage = 'Zero "Cons" parameter not allowed in simulation.'
            CALL PassMsg
            GOTO 898
         ENDIF
      ELSE                           ! in a "calibration" run
         IF (Cons .EQ. 0) THEN       ! allowed in a calibration run; recompute "Cons"
            AvgRR = 0.
            AvgHPlE = 0.
            DO I = 1, iNoDays
               Today = iFDMR + I - 1                ! todays date (sequence number)
               Tmin = Tmn(I) / 100.                 ! convert to degrees Celcius
               Tmax = Tmx(I) / 100.                 ! convert to degrees Celcius
               Precip = Prc(I) / 100000. * Area     ! convert to cubic meters
               DO J = 1, NumBC
                  IF (iBCD(J) .EQ. Today) THEN
                     IF (Snow(J) .GE. 0.) Snw  = Snow(J)
                  END IF
               END DO
               CALL SequenceDate(iD, iM, iY, Today); IF (ErrorLevel .NE. 0) GOTO 899
!
!              Compute daily insolation from minimum and maximum air temperatures and
!              from average mid-monthly cloudless daily solar insolation at the
!              surface.
!
               IF (iD .LE. 15) THEN
                  II = iM - 1
                  J = iY
                  IF (II .LT. 1) THEN
                     II = II + 12
                     J = J - 1
                  ENDIF
                  CALL DateSequence(15, II,  J,  J); IF (ErrorLevel .NE. 0) GOTO 899
                  CALL DateSequence(15, iM, iY, II); IF (ErrorLevel .NE. 0) GOTO 899
                  RR = (R(iM + 1) - R(iM)) / FLOAT(II - J) * FLOAT(Today - J) + R(iM)
               ELSE
                  II = iM + 1
                  J = iY
                  IF (II .GT. 12) THEN
                     II = II - 12
                     J = J + 1
                  END IF
                  CALL DateSequence(15, II, J, II); IF (ErrorLevel .NE. 0) GOTO 899
                  CALL DateSequence(15, iM, iY, J); IF (ErrorLevel .NE. 0) GOTO 899
                  RR = (R(iM + 2) - R(iM + 1)) / FLOAT(II - J) * FLOAT(Today - J) + R(iM + 1)
               END IF
               RR = RR * (B1 + B2 * AMIN1((TMAX - TMIN) / 15., 1.0)) * 10000. * AREA
!
!              Heat balance
!
               CALL SnowPackHeatBalance(AlbedS, Tmin, Tmax, Precip, Ta,   &
     &            Snw, Melt, Runoff)
               AvgRR = AvgRR + RR - Melt * 1000000. * 79.7
               AvgHPlE = AvgHPlE + HeatIndex(Ta, Tbase)
            END DO
            Cons = AvgRR / AvgHPlE   ! calories
            IF (Cons .LT. DPSILON) THEN
               Cons = 0
            ELSE
               CALL Round(Cons, 6)    ! this routine cannot handle a "zero"
            END IF
         END IF
      END IF
      cCons = Cons / Area / 10000.    ! always return the "Cons" value in cal/cm2.
!
!     Initialize volume summation variables.
!
      LTnsg = 0.
      LTrun = 0.
      LTint = 0.
      LTgwf = 0.
      LTuev = 0.
      LTlev = 0.
!
!     Daily loop
!
      DO I = 1, iNoDays
         Today = iFDMR + I - 1      ! todays date (sequence number)
!
!        Retrieve the daily met data for this day from the arrays.
!        Unit conversions: 
!           Temp from hundredths of degrees C -> degrees C
!           Precip from hundredths of mm -> cubic meters.
!
         Tmin = Tmn(I) / 100.
         Tmax = Tmx(I) / 100.
         Precip = Prc(I) / 100000. * Area
!
!        Use input boundary conditions, if available for this date.
!
         DO J = 1, NumBC
            IF (iBCD(J) .EQ. Today) THEN         
               IF (Upper(J)   .GE. 0.) Uszm = Upper(J)
               IF (Lower(J)   .GE. 0.) Lszm = Lower(J)
               IF (Ground(J)  .GE. 0.) Gzm  = Ground(J)
               IF (Surface(J) .GE. 0.) Ss   = Surface(J)
               IF (Snow(J)    .GE. 0.) Snw  = Snow(J)
            END IF
         END DO
!
!        Heat balance
!
         CALL SnowPackHeatBalance(AlbedS, Tmin, Tmax, Precip, Ta, Snw,    &
     &      Melt, Runoff)
         HPLSE = HeatIndex(Ta, Tbase) / (596. - .52 * Ta) / 1000000. * Cons
!
!        Mass balance  
!           (Runoff is in cubic meters when Outflow() finishes)
!
         LTnsg = LTnsg + Runoff  ! accumulate net supply to the ground surface
         CALL Outflow(Runoff) ; IF (ErrorLevel .NE. 0) GOTO 899
!
!        Save output.
!
         OutRunoff(I) = Runoff       ! cubic meters
         OutUszm(I)   = Uszm
         OutLszm(I)   = Lszm
         OutGzm(I)    = Gzm
         OutSs(I)     = Ss
         OutSnw(I)    = Snw
!
!        Accumulate flow volumes.
!
         LTrun = LTrun + vRun
         LTint = LTint + vInt
         LTgwf = LTgwf + vGW
         LTuev = LTuev + vUEv
         LTlev = LTlev + vLEv
      END DO
      RETURN
  898 ErrorLevel = 1
  899 RETURN
      END SUBROUTINE LBRM_Engine
!--------------------------------------------------------------------------
      REAL FUNCTION HeatIndex(Ta, Tbase)
      IMPLICIT NONE
      REAL, INTENT(IN) :: Ta, Tbase
      HeatIndex = EXP(Ta/Tbase)
!      IF (Ta .GE. 0) THEN
!         HeatIndex = Ta
!      ELSE
!         HeatIndex = 0.
!      ENDIF
      RETURN
      END FUNCTION HeatIndex
!--------------------------------------------------------------------------
      SUBROUTINE SnowPackHeatBalance(AlbedS, Tmin, Tmax, Precip, Ta, Snw, Melt, Runoff)
!
!     Snow pack heat balance.
!
      IMPLICIT NONE
      REAL, INTENT(IN) :: AlbedS, Tmin, Tmax, Precip
      REAL, INTENT(INOUT) :: Ta, Snw, Melt, Runoff
      REAL :: DD

      Ta = (Tmax + Tmin) / 2.
      Melt = 0.
      IF (Ta .GT. 0.) THEN
         IF (Snw .GE. 1.) THEN
            IF (Tmin .LT. 0.) THEN
              DD = Tmax ** 2 / (Tmax - Tmin) / 2.
            ELSE
              DD = Ta
            ENDIF
            Melt = AlbedS * DD
            IF (Melt .GT. Snw) Melt = Snw
         ENDIF
         Snw = Snw - Melt
         Runoff = Precip + Melt
      ELSE
         Snw = Snw + Precip
         Runoff = 0.
      ENDIF
      RETURN
      END SUBROUTINE SnowPackHeatBalance
!--------------------------------------------------------------------------
END MODULE LBRM_Main
!=============================================================================
