      subroutine irr_source(flag, j)

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine selects the irrigation source and performs the
!!    subtraction of irrigation water from source
!!
!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name           |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    aird(:)        |mm H2O        |amount of irrigation water applied to
!!                                  |HRU on current day that sums up to sol_sw
!!    auto_wstr(:)   |none or mm    |water stress threshold that triggers irrigation
!!    deepst(:)      |mm H2O        |depth of water in deep aquifer
!!    flowmin(:)     |m^3/s         |minimum instream flow for irrigation
!!                                  |diversions when IRR=1, irrigation water
!!                                  |will be diverted only when streamflow is
!!                                  |at or above FLOWMIN.
!!    hru_sub(:)     |none          |subbasin in which HRU is located
!!    inum1          |none          |reach or reservoir number
!!    irramt(:)      |mm H20        |depth of irrigation water applied to HRU
!!    irrno(:)       |none          |irrigation source location
!!                                  |if IRR=1, IRRNO is the number of the
!!                                  |          reach
!!                                  |if IRR=2, IRRNO is the number of the
!!                                  |          reservoir
!!                                  |if IRR=3, IRRNO is the number of the
!!                                  |          subbasin
!!                                  |if IRR=4, IRRNO is the number of the
!!                                  |          subbasin
!!                                  |if IRR=5, not used
!!    irrsc(:)       |none          |irrigation source code:
!!                                  |1 divert water from reach
!!                                  |2 divert water from reservoir
!!                                  |3 divert water from shallow aquifer
!!                                  |4 divert water from deep aquifer
!!                                  |5 divert water from unlimited source
!!    irr_mx(:)      |mm H2O        |maximum irrigation amount per auto application
!!    nhru           |none          |number of HRUs in watershed
!!    nirr(:)        |none          |sequence number of irrigation application
!!                                  |within the year
!!    res_vol(:)     |m^3           |reservoir volume
!!    rchstor(:)     |m^3           |water stored in reach
!!    rtwtr          |m^3           |water leaving reach on day
!!    shallst(:)     |mm H2O        |depth of water in shallow aquifer
!!    sirr           |string        |Stores the string for scheduled or auto-irrigation
!!    sol_sumfc(:)   |mm H2O        |amount of water held in the soil profile
!!                                  |at field capacity
!!    sol_sw(:)      |mm H2O        |amount of water stored in soil profile on any
!!                                  |given day
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!
!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    deepirr(:)  |mm H2O        |amount of water removed from deep aquifer
!!                               |for irrigation
!!    deepst(:)   |mm H2O        |depth of water in deep aquifer
!!    nirr(:)     |none          |sequence number of irrigation application
!!                               |within the year
!!    rchstor(:)  |m^3           |water stored in reach
!!    res_vol(:)  |m^3           |reservoir volume
!!    rtwtr       |m^3           |water leaving reach on day
!!    shallirr(:) |mm H2O        |amount of water removed from shallow aquifer
!!                               |for irrigation
!!    shallst(:)  |mm H2O        |depth of water in shallow aquifer
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!
!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    cnv         |none          |conversion factor (mm/ha => m^3)
!!    flag        |none          |irrigation flag:
!!                               |0 no irrigation operation on current day
!!                               |1 scheduled irrigation
!!                               |2 auto irrigation
!!    j           |none          |HRU number
!!    jrch        |none          |reach number
!!    jres        |none          |reservoir number
!!    k           |none          |counter
!!    vmm         |mm H2O        |amount of irrigation water applied to HRU
!!    vmma        |mm H2O        |maximum amount of water available in source
!!    vol         |m^3           |volume of irrigation water applied to HRU
!!    wtr_avail   |m^3           |volume of water available from reach
!!    wtrin       |m^3           |water outflow from reach prior to subtracting
!!                               |irrigation diversions
!!    xx          |???           |temporary variable to store temporary results
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!
!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    Intrinsic: Abs, Max, Min
!!    SWAT: irrigate
!!
!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~

      use parm
      implicit none

      integer :: flag, ii, j, jrch, jres, k
      logical :: filecheck
      real :: cnv, vmma, vmm, vol, wtr_avail, wtrin, xx

!! Clear variables
      cnv = 0.
      vmma = 0.
      vmm = 0.
      vol = 0.
      wtr_avail = 0.
      xx = 0.

!! Get identifier numbers
      jrch = 0
      jres = 0
      if (irrsc(j) == 1) jrch = inum1  !! If source = reach, get reach number
      if (irrsc(j) == 2) jres = inum1  !! If source = reservoir, get reservoir number

!! Fill the string for scheduled or auto-irrigation
      if (flag == 1) sirr(j) = "SCHD_IRR"
      if (flag == 2) sirr(j) = "AUTO_IRR"

     !! Calculate maximum irrigation amount =======================================================
      select case (irrsc(j))    !! Select irrigation source
        case (1)    !! Source --> Reach ---------------------------------------------------------
          if (irrno(j) == jrch) then  !! Check if water is to be removed from current reach
		        if (rtwtr > (flowmin(j) * 86400.)) then  !! Irrigate only if actual flow > minimum flow
		          cnv = 0.
		          cnv = hru_ha(j) * 10.
		          wtr_avail = 0.
		          wtr_avail = rtwtr + rchstor(jrch)
		          vmma = 0.  !! Calculate the maximum water available in source
		          vmma = (wtr_avail - flowmin(j) * 86400.) * flowfr(j) / cnv
		          if (divmax(j) < 0.) then  !! If divmax unit is 10^4 m^3, convert to mm H2O
		            xx = (Abs(divmax(j)) * 10000.) / cnv
		          else                      !! If divmax unit is mm H2O, get divmax
		            xx = divmax(j)
		          endif
		          vmma = Min(xx, vmma)      !! Calculate maximum irrigation water to be applied
		          xx = 0.  !! Clear temporary variable
		        end if
		      end if
        !! --------------------------------------------------------------------------------------

        case (2)    !! Source --> Reservoir -----------------------------------------------------
          if (irrno(j) == jres) then  !! Check if water is to be removed from current reservoir
		        cnv = 0.
		        cnv = hru_ha(j) * 10.
		        vmma = 0.
		        vmma = res_vol(jres) / cnv     !! Calculate the maximum water available in source
          end if
        !! --------------------------------------------------------------------------------------

        case (3)    !! Source --> Shallow aquifer -----------------------------------------------
          vmma = 0.
          do k = 1, nhru    !! Loop throughout HRUs to get water available from shallow aquifer
            if (hru_sub(k) == irrno(j)) then
              cnv = 0.
              cnv = hru_ha(k) * 10.
              vmma = vmma + shallst(k) * cnv  !! Calculate the maximum water available in source
            end if
          end do
          cnv = 0.
          cnv = hru_ha(j) * 10.
          vmma = vmma / cnv              !! Calculate maximum irrigation water to be applied
        !! --------------------------------------------------------------------------------------

        case (4)   !! Source --> Deep aquifer ---------------------------------------------------
          vmma = 0.
          do k = 1, nhru    !! Loop throughout HRUs to get the water available from deep aquifer
            if (hru_sub(k) == irrno(j)) then
              cnv = 0.
              cnv = hru_ha(k) * 10.
              vmma = vmma + deepst(k) * cnv  !! Calculate the maximum water available in source
            end if
          end do
          cnv = 0.
          cnv = hru_ha(j) * 10.
          vmma = vmma / cnv              !! Calculate maximum irrigation water to be applied
        !! --------------------------------------------------------------------------------------

         case (5)   !! Source - Unlimited --------------------------------------------------------
          vmma = 0.
          vmma = sol_sumfc(j)             !! Calculate maximum irrigation water to be applied
         !! --------------------------------------------------------------------------------------

      end select !! =============================================================================



!! Check if there is water available from source and, if so, apply irrigation operation
      if (vmma > 0.) then
        vmm = 0.
        vmm = Min(sol_sumfc(j), vmma)  !! Check if available water from source is
                                       !!   lower than amount of water held in the soil
                                       !!   profile at field capacity
     !! Check if available water from source < maximum irrigation amount
        if (flag == 1) vmm = Min(vmm, irramt(j))  !! Scheduled irrigation
        if (flag == 2) vmm = Min(vmm, irr_mx(j))  !! Auto-irrigation

     !! Subtract irrigation water applied from source ===========================================
        select case (irrsc(j))     !! Select irrigation source
          case (1)    !! Source --> Reach -------------------------------------------------------
            if (irrno(j) == jrch) then  !! Check if water is to be removed from current reach
		          wtrin = 0.
		          wtrin = rtwtr + rchstor(jrch)  !! Update water outflow from reach before irrigation diversions
		          cnv = 0.
		          cnv = hru_ha(j) * 10.
		          vol = 0.
		          vol = vmm * cnv                !! Irrigation amount in m^3/d
		          vol = Min(vol, wtr_avail)      !! Check if irrigation amount > water available
		          if (rchstor(jrch) > vol) then  !! Remove water first from rchstor
		            xx = 0.
		            rchstor(jrch) = rchstor(jrch) - vol
		          else                           !! If water at rchstor is not enough, remove the
		            xx = vol - rchstor(jrch)     !!   difference from rtwrt
		            vol = rchstor(jrch)
		            rchstor(jrch) = 0.
		            rtwtr = rtwtr - xx
		            if (rtwtr < 0.) then    !! If source volume gets below zero, fix it
		            	xx = xx + rtwtr
		            	rtwtr = 0.
		            end if
		         !! Subtract irrigation from reach outflow
		            if (ievent > 0) then
		              do ii = 1, nstep
		                hrtwtr(ii) = hrtwtr(ii) - xx
		                hrtwtr(ii) = Max(0., hrtwtr(ii))
		              end do
		            end if
		          end if
		          vol = vol + xx
		          vmm = vol / cnv
		          xx = 0.  !! Clear temporary variable
            end if
          !! ------------------------------------------------------------------------------------

          case (2)    !! Source --> Reservoir ---------------------------------------------------
            if (irrno(j) == jres) then  !! Check if water is to be removed from current reservoir
		          cnv = 0.
		          cnv = hru_ha(j) * 10.
		          vol = 0.
		          vol = vmm * cnv    !! Irrigation amount in m^3/d
		          res_vol(jres) = res_vol(jres) - vol
		          if (res_vol(jres) < 0.) then    !! If source volume gets below zero, fix it
		            vol = vol + res_vol(jres)
		            res_vol(jres) = 0.
		          end if
		          vmm = vol / cnv
            end if
          !! ------------------------------------------------------------------------------------

          case (3)    !! Source --> Shallow aquifer ---------------------------------------------
            do k = 1, nhru    !! Loop throughout HRUs to get water available from shallow aquifer
              if (hru_sub(k) == irrno(j)) then
                shallst(k) = shallst(k) - vmm
                if (shallst(k) < 0.) then    !! If source volume gets below zero, fix it
                  vmm = vmm + shallst(k)
                  shallst(k) = 0.
                end if
                shallirr(k) = shallirr(k) + vmm
              end if
            end do
          !! ------------------------------------------------------------------------------------

          case (4)   !! Source --> Deep aquifer -------------------------------------------------
            do k = 1, nhru    !! Loop throughout HRUs to get water available from deep aquifer
              if (hru_sub(k) == irrno(j)) then    !! If source volume gets below zero, fix it
                deepst(k) = deepst(k) - vmm
                if (deepst(k) < 0.) then
                  vmm = vmm + deepst(k)
                  deepst(k) = 0.
                end if
                deepirr(k) = deepirr(k) + vmm
              end if
            end do
          !! ------------------------------------------------------------------------------------

        end select !! ===========================================================================

     !! Calculate the irrigation amount in m^3/d
        cnv = 0.
        cnv = hru_ha(j) * 10.
        vol = 0.
        vol = vmm * cnv    !! Irrigation amount in m^3/d

     !! Calculate the irrigation amount and losses
        call irrigate(flag, j, vmm)

     !! Check if HRU is a pothole
        if (pot_fr(j) > 1.e-6) then
          pot_vol(j) = pot_vol(j) + aird(j) / (10. * potsa(j))
        end if

     !! Advance irrigation operation number
        nirr(j) = nirr(j) + 1

     !! Write the outputs to 'output.mgt' file
        if (imgt == 1) then
          write (143, 1000) subnum(j), hruno(j), iyr, i_mo, iida,
     *    hru_km(j), "         ",  sirr(j), phubase(j), phuacc(j),
     *    sol_sw(j), bio_ms(j), sol_rsd(1,j),sol_sumno3(j),
     *    sol_sumsolp(j), aird(j), irrsc(j), irrno(j)
        end if

     !! Write the outputs to 'output.irr' file
      inquire(file='output.irr', exist=filecheck)
      if (filecheck) then
        open(9143, file='output.irr', action='write', status='old',
     *       access='append')
      else
        open(9143, file='output.irr', action='write', status='replace')
        write(9143, 1143)
      end if
      write (9143, '(a5)', advance = 'no') subnum(j)
      write (9143, '(a1)', advance = 'no') ","
      write (9143, '(a4)', advance = 'no') hruno(j)
      write (9143, '(a1)', advance = 'no') ","
      write (9143, '(i6)', advance = 'no') iyr
      write (9143, '(a1)', advance = 'no') ","
      write (9143, '(i6)', advance = 'no') i_mo
      write (9143, '(a1)', advance = 'no') ","
      write (9143, '(i6)', advance = 'no') iida
      write (9143, '(a1)', advance = 'no') ","
      write (9143, '(e10.5)', advance = 'no') hru_km(j)
      write (9143, '(a1)', advance = 'no') ","
      write (9143, '(a15)', advance = 'no') sirr(j)
      write (9143, '(a1)', advance = 'no') ","
      write (9143, '(f10.2)', advance = 'no') phubase(j)
      write (9143, '(a1)', advance = 'no') ","
      write (9143, '(f10.2)', advance = 'no') sol_sw(j)
      write (9143, '(a1)', advance = 'no') ","
      write (9143, '(f10.2)', advance = 'no') tird(j)
      write (9143, '(a1)', advance = 'no') ","
      write (9143, '(f10.2)', advance = 'no') aird(j)
      write (9143, '(a1)', advance = 'no') ","
      write (9143, '(f10.2)', advance = 'no') qird(j)
      write (9143, '(a1)', advance = 'no') ","
      write (9143, '(f10.2)', advance = 'no') lird(j)
      write (9143, '(a1)', advance = 'no') ","
      write (9143, '(i10)', advance = 'no') irrsc(j)
      write (9143, '(a1)', advance = 'no') ","
      write (9143, '(i10)', advance = 'no') irrno(j)
      write (9143, *)
      close (9143)

     !! If irrigation source is reach, update reach specific parameters
		    if (irrsc(j) == 1) then
		      if (wtrin /= rtwtr .and. wtrin > 0.01) then
		        sedrch  = sedrch  * rtwtr / wtrin
		        rch_san = rch_san * rtwtr / wtrin
		        rch_sil = rch_sil * rtwtr / wtrin
		        rch_cla = rch_cla * rtwtr / wtrin
		        rch_sag = rch_sag * rtwtr / wtrin
		        rch_lag = rch_lag * rtwtr / wtrin
		        rch_gra = rch_gra * rtwtr / wtrin
		        if (sedrch  < 1.e-6) then
		          sedrch  = 0.
		          rch_san = 0.
		          rch_sil = 0.
		          rch_cla = 0.
		          rch_sag = 0.
		          rch_lag = 0.
		          rch_gra = 0.
		        end if
		        if (ievent > 0) then
		          do ii = 1, nstep
		            hsedyld(ii) = hsedyld(ii) * rtwtr / wtrin
		          end do
		        end if
		      end if
		    end if

      end if

1000  format (a5,1x,a4,3i6,1x,e10.5,1x,2a15,7f10.2,10x,f10.2,70x,i10,
     *        10x,i10)
1143  format ("SUB,HRU,YEAR,MONTH,DAY,HRU_AREA,IRR,PHUbase,SW,"//
     +        "IRRt,IRRf,IRRq,IRRs,IRRsc,IRRno")

      return
      end
