C
C*************************** IRR_READ_ONE ******************************
C
C Read in irrigation data for the first time step (times 0 and DELTAT). 
C
C Description of the irrigation management scheme:
C
C Irrigation rates (in units [L/T]) are read in from input file unit IINDEM12,
C together with various control parameters:
C IRR_FLAG (integer scalar)
C   = 1 to activate the irrigation module.
C   = 0 to skip all irrigation input and processing.
C IRR_SW (real scalar, units [-])
C   The threshold value of water saturation (SW = moisture content/porosity)
C   that determines whether available irrigation water is applied or not. I.e.,
C   at any surface node I during the course of a simulation, any available
C   irrigation water at node I is applied (i.e., is added to ATMPOT(I)) only if
C   SW(I) is less than IRR_SW.
C   This very rudimentary scheme for "irrigation management" can be readily
C   extended to more complex schemes in future work (space and time varying
C   threshold SW values, threshold values that are dependent on soil properties,
C   an incremental rather than threshold criterion for applying irrigation
C   water, etc).
C IRR_HT (integer scalar)
C   = 1 for constant in time irrigation rates.
C   = 0 for temporally variable irrigation rates.
C   IRR_HT is not actually read in, but rather it is set in this subroutine and
C   in subroutine IRR_READ_NXT according to whether there are (still) TIMEIN
C   values to be read in from unit IINDEM12. I.e., as soon as there are no
C   longer TIMEIN values in the input file, IRR_HT is set to 1 and the
C   irrigation rates remain constant in time (and equal to those read in at the
C   previous TIMEIN value) for the rest of the simulation.
C IRR_HS (integer scalar)
C   = 1 for spatially homogeneous irrigation rates.
C   = 0 for spatially variable irrigation rates.
C   In the case IRR_HS = 0, a raster map of irrigation rates is read in. In the
C   case IRR_HS = 1, a single irrigation rate is read in and is copied to the
C   entire raster.
C
C The input raster map of irrigation rates follows the matrix convention, just
C as is done with the DEM_MAP, ZONE, BASE_MAP, LAKES_MAP, and VEG_MAP inputs
C (and in contrast to all the "DTM_xxx" input raster maps that are generated by
C the preprocessor and that are read in according to a Cartesian reference
C frame). 
C
C Note that in this implementation, the irrigation module can only be used in
C "DEM input" mode (i.e., "GRID input" mode is not supported) and with the
C subsurface flow module activated.
C
C Other variables specific to the irrigation module:
C IRR_TIM (real array, dimension 3, units [T])
C   Most current input TIMEIN values for the irrigation data, with
C   IRR_TIM(1) < IRR_TIM(2) < IRR_TIM(3) and IRR_TIM(2) < TIME <= IRR_TIM(3),
C   where TIME is the time at the current simulation time level.
C   The reason that three consecutive TIMEIN values need to be stored instead
C   of just two is to properly handle back-stepping occurrences.
C   [Note: back-step handling for the irrigation scheme still needs to be
C   thoroughly tested, as it should probably be done a bit differently than
C   how it is handled in the case of atmospheric BC input and non-atmospheric
C   Neumann and Dirichlet BC input. In the case of these BCs, all processing is
C   done at the end of a time step, whereas for irrigation data the processing
C   is done at the beginning of a time step.]
C   During the simulation, at each TIME value (where as just mentioned we have
C   IRR_TIM(2) < TIME <= IRR_TIM(3)), the irrigation rate for each cell in the
C   raster map is the value that was read in at TIMEIN = IRR_TIM(2). In other
C   words, the input irrigation data is interpreted as a piecewise constant
C   function, or ietograph.
C IRR_INP_RAST (real array, dimension 3 x ROWMAX x COLMAX, units [L/T])
C   The raster map of input irrigation rates (see also the description of IRR_HS
C   above) at the three IRR_TIM time levels.
C IRR_NELT (real array, dimension NTRMAX, units [L/T])
C   The irrigation rates at the current TIME level, as transferred from the
C   cells of the surface DEM (raster map input) to the triangles of the surface
C   grid for the subsurface module. In the case of DOSTEP = 1, the irrigation
C   rate for each cell is simply assigned to each of the two triangles that are
C   formed from that cell. In the case of DOSTEP > 1 (grid coarsening), the
C   transfer is rather more elaborate and also involves a (simple) averaging
C   procedure. See subroutine IRR_DEMCELL_TO_NODE.
C IRR_NNOD (real array, dimension NODMAX, units [L/T])
C   Ultimately, the irrigation rates are needed at the surface nodes since it
C   is at this level that atmospheric inputs are processed (rainfall-runoff
C   partitioning, boundary condition switching, etc). This transfer from
C   elements (triangles) to nodes, again involving simple averaging, is also
C   done in subroutine IRR_DEMCELL_TO_NODE.
C   Note that there is actually one final transfer needed: from [L/T] units to
C   [L^3/T]. This is not done until available irrigation water is added to
C   ATMPOT (only for nodes that don't exceed the IRR_SW threshold). This final
C   step is performed in subroutines ATMONE, ATMNXT, and ATMBAK.
C IRR_ATMPOT_SAV (real array, dimension NODMAX, units [L^3/T])
C   ATMPOT values need to be saved before any irrigation water is added, so
C   that at the next time step we don't add irrigation water to an ATMPOT
C   array that has already had irrigation water added to it. This needs to
C   be done because the input TIMEIN values for atmospheric forcing BCs will
C   not generally coincide with the TIMEIN values for irrigation data input.
C IRR_POT_ATM, IRR_POT_IRRWATER,
C IRR_POT_ATM_IRR, IRR_IRRWATER_NOT_USED (real scalars, units [L^3/T])
C   Variables for IOUT33 output (see below).
C IRR_POT_ATM_CUMV, IRR_POT_IRRWATER_CUMV,
C IRR_POT_ATM_IRR_CUMV, IRR_IRRWATER_NOT_USED_CUMV (real scalars, units [L^3])
C   Additional variables for IOUT33 output (see below).
C
C All variables specific to the irrigation module are defined in the IRR_CBLOCK
C common block.
C
C The output file unit for the irrigation scheme is IOUT33 and consists of a
C hydrograph-type output of irrigation rate and cumulative volume and how this
C irrigation water is partitioned. Additional irrigation-related tracking
C (e.g., vtk-type output) can be added as needed.
C
C***********************************************************************
C
      SUBROUTINE IRR_READ_ONE(NROW,NCOL,NNOD,NTRI,DOSTEP,
     1                        TIME,DELTAT,
     2                        TRIANG,TP2D,DEM_MAP)
C
      IMPLICIT  NONE
      INCLUDE  'CATHY.H'
      INCLUDE  'IOUNITS.H'
      INCLUDE  'IRR_CBLOCK.H'
      INTEGER   I,J,IRR_I
      INTEGER   NROW,NCOL,NNOD,NTRI,DOSTEP
      INTEGER   TRIANG(4,*),TP2D(*)
      REAL*8    TIMEIN,KRATE
      REAL*8    TIME,DELTAT
      REAL*8    DEM_MAP(ROWMAX,*)
C     
      IRR_HT = 0
      IRR_TIM(1) = 0.0D0
      IRR_TIM(2) = 0.0D0
      DO I = 1,NNOD
         IRR_NNOD(I) = 0.0D0
      END DO
      DO I = 1,NROW
         DO J = 1,NCOL
            IRR_INP_RAST(1,I,J) = 0.0D0
            IRR_INP_RAST(2,I,J) = 0.0D0
         END DO
      END DO
      READ(IINDEM12,*) IRR_SW
      READ(IINDEM12,*) IRR_HS
      READ(IINDEM12,*,END=600) TIMEIN
      IF (DELTAT .GE. 1.0D+15) TIMEIN = 0.0D0
      IRR_TIM(3)=TIMEIN
      IF (IRR_HS .EQ. 0) THEN
         DO I = 1,NROW
            READ(IINDEM12,*) (IRR_INP_RAST(3,I,J),J=1,NCOL)
         END DO
      ELSE
         IRR_HS = 1
         READ(IINDEM12,*) KRATE
         DO I = 1,NROW
            DO J = 1,NCOL
               IRR_INP_RAST(3,I,J) = KRATE
            END DO
         END DO
      END IF
      IF (DELTAT .GE. 1.0D+15) GO TO 300
  200 IF (TIME .LE. IRR_TIM(3)) GO TO 300
      IRR_TIM(1)=IRR_TIM(2)
      IRR_TIM(2)=IRR_TIM(3)
      DO I = 1,NROW
         DO J = 1,NCOL
            IRR_INP_RAST(1,I,J) = IRR_INP_RAST(2,I,J)
            IRR_INP_RAST(2,I,J) = IRR_INP_RAST(3,I,J)
         END DO
      END DO
      READ(IINDEM12,*,END=400) TIMEIN
      IRR_TIM(3)=TIMEIN
      IF (IRR_HS .EQ. 0) THEN
         DO I = 1,NROW
            READ(IINDEM12,*) (IRR_INP_RAST(3,I,J),J=1,NCOL)
         END DO
      ELSE
         READ(IINDEM12,*) KRATE
         DO I = 1,NROW
            DO J = 1,NCOL
               IRR_INP_RAST(3,I,J) = KRATE
            END DO
         END DO
      END IF
      GO TO 200
  300 CONTINUE
      IF (IRR_TIM(3) .GT. IRR_TIM(2)) THEN
         IRR_I = 2
      ELSE
         IRR_I = 3
      END IF
      CALL IRR_DEMCELL_TO_NODE(NROW,NCOL,NNOD,NTRI,DOSTEP,IRR_I,
     1                         TRIANG,TP2D,DEM_MAP)
      GO TO 800
C
  400 IRR_HT = 1
      GO TO 300
C
C If there is no input value at all of TIMEIN, then this is equivalent to
C irrigation rates being zero for the entire simulation (constant in time), 
C in which case it is more efficient to simply avoid all irrigation processing,
C so we set IRR_FLAG to zero.
C
  600 IRR_FLAG = 0
C
  800 RETURN
      END
