!! manage water intakes from river for irrigation purposes
!|author:  <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a>
! license: <a href="http://www.gnu.org/licenses/">GPL</a>
!    
!### History
!   
! current version  1.9 - 19th March 2025  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 27/Jun/2012 | Original code |
! | 1.1      | 10/Oct/2023 | an irrigated area can accept multiple sources from different intakes|
! | 1.2      | 13/Nov/2023 | changed District to Intake and eta (irrigation efficiency) added |
! | 1.3      | 04/Dec/2023 | maximum discharge and eflow are defined daily. function SetDailyArray added |
! | 1.4      | 06/May/2024 | id written as header in output file instead of name |
! | 1.5      | 04/Jun/2024 | output file for both diverted and downstream discharge |
! | 1.6      | 05/Jun/2024 | multiple intakes in one cell are allowed |
! | 1.7      | 31/Jul/2024 | output files are not created when dtOut = 0 |
! | 1.8      | 21/Aug/2024 | upstream discharge written to output file |
! | 1.9      | 19/Mar/2025 | fluxQ set to zero when out of irrigation period |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! routines to model water intakes from river for irrigation purposes
!
MODULE Irrigation         
			
! Modules used: 

USE DataTypeSizes, ONLY: &
    !Imported parameetrs:
    short, float

USE GeoLib, ONLY: &
    !Imported definitions:
    Coordinate

USE GridLib, ONLY: &
    !Imported definitions:
    grid_integer, grid_real, &
    !Imported routines:
    NewGrid


USE GridOperations, ONLY: &
    !Imported routines:
    GridByIni, GetArea, &
    GetIJ, ASSIGNMENT(=)

USE GridStatistics, ONLY: &
    !imported routines:
    GetMean

USE IniLib, ONLY: &
    !Imported definitions:
    IniList, &
    !Imported routines:
    IniOpen, IniReadInt, IniReadString, &
    IniReadReal, KeyIsPresent

USE StringManipulation, ONLY: &
    !Imported routines:
    ToString, StringToFloat

USE LogLib, ONLY: &
    !Imported routines:
    Catch

USE Chronos, ONLY: &
    !Imported definitions:
    DateTime, timeString, &
    !Imported routines:
    DayOfYear, &
    !imported operators:
    ASSIGNMENT (=), &
    OPERATOR (+), &
    OPERATOR (==)
    

USE Utilities, ONLY: &
    !imported routines:
    GetUnit

USE TableLib, ONLY : &
    !Imported definitions:
    Table, &
    !Imported routines:
    TableNew, TableGetNrows, TableGetValue

USE DomainProperties, ONLY : &
    !Imported variables:
    mask

IMPLICIT NONE

!Type definitions:

TYPE Intake
  CHARACTER (LEN = 300)  :: name
  CHARACTER (LEN = 300)  :: id
  REAL (KIND = float)    :: max_discharge (365) !! maximum concessed discharge for every day of the year (m3/s)
  REAL (KIND = float)    :: fluxQ !! actual diverted discharge (m3/s)
  TYPE (Coordinate)      :: xy !!easting and northing 
  REAL (KIND = float)    :: e_flow (365) !! minimum environmental flow for every day of the year (m3/s)
  REAL (KIND = float)    :: eta !! irrigation efficiency (0-1)
  INTEGER (KIND = short) :: doy_start !!day of year intake starts taking water
  INTEGER (KIND = short) :: doy_stop !!day of year intake stops taking water
  INTEGER (KIND = short) :: r, c !!row and column number in the local coordinate of raster domain
  REAL (KIND = float)    :: sat_max !!district averaged  soil saturation above which irrigation is stopped
  TYPE (grid_integer)    :: mask
  REAL (KIND = float)    :: area !! surface (m2)
END TYPE Intake
 

TYPE Intakes
   INTEGER (KIND = short) :: count !!number of intakes
   TYPE (Intake), ALLOCATABLE :: elem (:)
END TYPE Intakes

PUBLIC :: IrrigationConfig
PUBLIC :: IrrigationUpdate

TYPE (grid_real)       :: Qirrigation
TYPE (grid_real)       :: irrigationFlux
INTEGER (KIND = short) :: dtIrrigation
TYPE (Intakes)         :: fields


INTEGER (KIND = short), PRIVATE :: unitIrrigationDownstream
INTEGER (KIND = short), PRIVATE :: unitIrrigationUpstream
INTEGER (KIND = short), PRIVATE :: unitIrrigationDiverted 
CHARACTER (LEN = 15), PRIVATE :: epsg

!private declarations

TYPE (grid_real), PRIVATE :: cpQriver !! local copy of river discharge



!PRIVATE routines:
PRIVATE :: SetDailyArray

!=======         
CONTAINS
!======= 
! Define procedures contained in this module. 

!==============================================================================
!| Description:
!   configure irrigation
SUBROUTINE IrrigationConfig &
!
(file, path_out, dtOut)

IMPLICIT NONE

!Arguments with intent (in)
CHARACTER (LEN = *), INTENT (IN) :: file
CHARACTER (LEN = *), INTENT(IN) :: path_out
INTEGER (KIND = short), INTENT(IN) :: dtOut !!dt out irrigation

!local declarations
TYPE (IniList) :: iniDb
CHARACTER (LEN = 10)   :: secNumber
INTEGER (KIND = short) :: k
CHARACTER (LEN = 300) :: string
LOGICAL :: check

!-----------------------------------end of declarations------------------------

CALL IniOpen (file, iniDB)

!read EPSG (used only for witing string into output file)
epsg = IniReadString ('epsg', iniDB)

!Allocate fields
fields % count = IniReadInt ('count', iniDB)
ALLOCATE ( fields % elem (fields % count) )

!read intake and irrigation district properties
DO k = 1, fields % count
  secNumber = ToString (k)
  fields % elem (k) % name = IniReadString ('name', iniDB, section = secNumber)
  fields % elem (k) % id = IniReadString ('id', iniDB, section = secNumber)
  string = IniReadString ('max-discharge', iniDB, section = secNumber)
  fields % elem (k) % max_discharge = SetDailyArray (string)
  fields % elem (k) % xy % northing = IniReadReal ('northing', &
                                      iniDB, section = secNumber)
  fields % elem (k) % xy % easting = IniReadReal ('easting', iniDB, &
                                     section = secNumber)
  string = IniReadString ('e-flow', iniDB, section = secNumber)
  fields % elem (k) % e_flow = SetDailyArray (string)
  fields % elem (k) % doy_start = IniReadInt ('doy-start', iniDB, &
                                  section = secNumber)
  fields % elem (k) % doy_stop = IniReadInt ('doy-stop', iniDB, &
                                 section = secNumber)
  fields % elem (k) % sat_max = IniReadReal ('sat-max', iniDB, &
                                section = secNumber)
  !read irrigation efficiency. Default = 1
  IF (KeyIsPresent ('eta', iniDB, section = secNumber)) THEN
      fields % elem (k) % eta = IniReadReal ('eta', iniDB, &
                                section = secNumber)
      IF ( fields % elem (k) % eta > 1. ) THEN
          CALL Catch ('warning', 'Irrigation',  &
          'irrigation efficiency > 1 in intake: ',  &
           argument = fields % elem (k) % name )
      END IF
      IF ( fields % elem (k) % eta < 0. ) THEN
          CALL Catch ('warning', 'Irrigation',  &
          'irrigation efficiency < 0 in intake: ',  &
           argument = fields % elem (k) % name )
      END IF
  ELSE
      fields % elem (k) % eta = 1.
  END IF
  !load i_th subbasin mask
  CALL GridByIni (iniDB, fields % elem (k) % mask, secNumber, 'mask')
  fields % elem (k) % area = GetArea (fields % elem (k) % mask)
  
  !set intake location in local grid reference system
  CALL GetIJ (fields % elem (k) % xy % easting, &
              fields % elem (k) % xy % northing, &
              mask, fields % elem (k) % r, fields % elem (k) % c, check)
        IF (.NOT. check) THEN
          CALL Catch ('error', 'Irrigation',  &
             'river intake located out of discharge grid in intake: ',  &
              argument = TRIM(fields % elem (k) % name) )
        END IF
END DO

!Open files for output

IF ( dtOut > 0 ) THEN

    !diverted discharge
    unitIrrigationDiverted = GetUnit ()
    OPEN (unit = unitIrrigationDiverted, &
          file = path_out(1:LEN_TRIM(path_out))//'irrigation_diverted.fts')

    !populate output file with metadata information
    WRITE (unitIrrigationDiverted,'(a)') 'description = irrigation discharge &
                                          diverted from water courses'
    WRITE (unitIrrigationDiverted,'(a)') 'unit = m3/s'
    WRITE (unitIrrigationDiverted,'(a)') 'epsg = ' // TRIM (epsg)
    WRITE (unitIrrigationDiverted,'(a,i)') 'count = ', fields % count
    WRITE (unitIrrigationDiverted,'(a,i)') 'dt = ', dtOut
    WRITE (unitIrrigationDiverted,'(a)') 'missing-data =   -999.900'
    WRITE (unitIrrigationDiverted,'(a)') 'offsetz =      0.000'

    !metadata section
    WRITE (unitIrrigationDiverted,*) !blank line
    WRITE (unitIrrigationDiverted,'(a)') 'metadata'

    DO k = 1, fields % count
        WRITE (unitIrrigationDiverted,'(a,3x,a,3x,3f15.3)') &
              TRIM(fields % elem (k) % name), &
              TRIM(fields % elem (k) % id), &
              fields % elem (k) % xy % easting, &
              fields % elem (k) % xy % northing, 0.00
    END DO

    !data section
    WRITE (unitIrrigationDiverted,*) !blank line
    WRITE (unitIrrigationDiverted,'(a)') 'data'
    WRITE (unitIrrigationDiverted,'(a)', advance = 'no') 'DateTime '

    DO k = 1, fields % count - 1
        WRITE (unitIrrigationDiverted,'(a, 1x)', advance = 'no') &
               TRIM (fields % elem (k) % id)
    END DO

    WRITE (unitIrrigationDiverted,'(a)') &
           TRIM ( fields % elem (fields % count) % id )


    !downstream discharge
    unitIrrigationDownstream = GetUnit ()
    OPEN (unit = unitIrrigationDownstream, &
          file = path_out(1:LEN_TRIM(path_out))//'irrigation_downstream.fts')

    !populate output file with metadata information
    WRITE (unitIrrigationDownstream,'(a)') 'description = irrigation  &
                                            discharge downstream the intake'
    WRITE (unitIrrigationDownstream,'(a)') 'unit = m3/s'
    WRITE (unitIrrigationDownstream,'(a)') 'epsg = ' // TRIM (epsg)
    WRITE (unitIrrigationDownstream,'(a,i)') 'count = ', fields % count
    WRITE (unitIrrigationDownstream,'(a,i)') 'dt = ', dtOut
    WRITE (unitIrrigationDownstream,'(a)') 'missing-data =   -999.900'
    WRITE (unitIrrigationDownstream,'(a)') 'offsetz =      0.000'

    !metadata section
    WRITE (unitIrrigationDownstream,*) !blank line
    WRITE (unitIrrigationDownstream,'(a)') 'metadata'

    DO k = 1, fields % count
        WRITE (unitIrrigationDownstream,'(a,3x,a,3x,3f15.3)') &
              TRIM(fields % elem (k) % name), &
              TRIM(fields % elem (k) % id), &
              fields % elem (k) % xy % easting, &
              fields % elem (k) % xy % northing, 0.00
    END DO

    !data section
    WRITE (unitIrrigationDownstream,*) !blank line
    WRITE (unitIrrigationDownstream,'(a)') 'data'
    WRITE (unitIrrigationDownstream,'(a)', advance = 'no') 'DateTime '

    DO k = 1, fields % count - 1
        WRITE (unitIrrigationDownstream,'(a, 1x)', advance = 'no') &
               TRIM (fields % elem (k) % id)
    END DO

    WRITE (unitIrrigationDownstream,'(a)') &
           TRIM ( fields % elem (fields % count) % id )
    
    
    !upstream discharge
    unitIrrigationUpstream = GetUnit ()
    OPEN (unit = unitIrrigationUpstream, &
          file = path_out(1:LEN_TRIM(path_out))//'irrigation_upstream.fts')

    !populate output file with metadata information
    WRITE (unitIrrigationUpstream,'(a)') 'description = irrigation  &
                                            discharge upstream the intake'
    WRITE (unitIrrigationUpstream,'(a)') 'unit = m3/s'
    WRITE (unitIrrigationUpstream,'(a)') 'epsg = ' // TRIM (epsg)
    WRITE (unitIrrigationUpstream,'(a,i)') 'count = ', fields % count
    WRITE (unitIrrigationUpstream,'(a,i)') 'dt = ', dtOut
    WRITE (unitIrrigationUpstream,'(a)') 'missing-data =   -999.900'
    WRITE (unitIrrigationUpstream,'(a)') 'offsetz =      0.000'

    !metadata section
    WRITE (unitIrrigationUpstream,*) !blank line
    WRITE (unitIrrigationUpstream,'(a)') 'metadata'

    DO k = 1, fields % count
        WRITE (unitIrrigationUpstream,'(a,3x,a,3x,3f15.3)') &
              TRIM(fields % elem (k) % name), &
              TRIM(fields % elem (k) % id), &
              fields % elem (k) % xy % easting, &
              fields % elem (k) % xy % northing, 0.00
    END DO

    !data section
    WRITE (unitIrrigationUpstream,*) !blank line
    WRITE (unitIrrigationUpstream,'(a)') 'data'
    WRITE (unitIrrigationUpstream,'(a)', advance = 'no') 'DateTime '

    DO k = 1, fields % count - 1
        WRITE (unitIrrigationUpstream,'(a, 1x)', advance = 'no') &
               TRIM (fields % elem (k) % id)
    END DO

    WRITE (unitIrrigationUpstream,'(a)') &
           TRIM ( fields % elem (fields % count) % id )

END IF

!allocate variables:
CALL NewGrid (Qirrigation, mask, 0.)
CALL NewGrid (irrigationFlux, mask, 0.)
CALL NewGrid (cpQriver, mask, 0.)


RETURN
END SUBROUTINE IrrigationConfig


!==============================================================================
!| Description:
!   update irrigation and write output
SUBROUTINE IrrigationUpdate &
!
(time, sat, Qriver, dtOut, timeOut)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime),  INTENT (IN) :: time !!current time
TYPE (grid_real), INTENT (IN) :: sat !!soil saturation (0-1)
TYPE (grid_real), INTENT (IN) :: Qriver  !!river discharge at current time
INTEGER (KIND = short), INTENT(IN) :: dtOut !!time step for output

!Argument with intent (inout):
TYPE (DateTime), INTENT(INOUT) :: timeOut

!Local declarations:
INTEGER (KIND = short) :: i,j,k,r,c
REAL (KIND = float) :: flux !!irrigation flux (m/s)
REAL (KIND = float) :: fluxQ !!irrigation disharge (m3/s)
INTEGER (KIND = short) :: doy
REAL (KIND = float)    :: meanSat
REAL (KIND = float) :: Qdownstream

!----------------------------end of declarations-------------------------------

timeString = time

irrigationFlux = 0.

Qirrigation = 0.

cpQriver = Qriver

DO k = 1, fields % count 
    !check soil saturation and irrigation season
    doy = DayOfYear (time)
    fluxQ = 0.
    flux = 0.
    !meanSat = GetMean (sat, maskInteger = fields % elem (k) % mask )
    IF (doy > fields % elem (k) % doy_start .AND. &
        doy < fields % elem (k) % doy_stop  ) THEN ! .AND. &
        !meanSat < fields % elem (k) % sat_max) THEN
        !compute irrigation flux [m/s]
        r = fields % elem (k) % r
        c = fields % elem (k) % c
        
 
        IF (cpQriver % mat (r,c) > fields % elem (k) % e_flow (doy) ) THEN
           fluxQ = cpQriver % mat (r,c) - fields % elem (k) % e_flow (doy)
           !limit to concessed discharge
           IF (fluxQ > fields % elem (k) % max_discharge (doy) ) THEN 
             fluxQ = fields % elem (k) % max_discharge (doy)
           END IF
           
           !store irrigation discharge
           Qirrigation % mat (r,c) =  Qirrigation % mat (r,c) + fluxQ
           
           !conversion m3/s -> m/s
           flux = fluxQ / fields % elem (k) % area  
           
        ELSE
           fluxQ = 0.
           flux = 0. 
        END IF
        
        fields % elem (k) % fluxQ = fluxQ
        
        cpQriver % mat (r,c) = cpQriver % mat (r,c) - fluxQ
        
        !populate irrigation map 
        DO i = 1, irrigationFlux % idim
          DO j = 1, irrigationFlux % jdim
            IF (fields % elem (k) % mask % mat(i,j) /= &
                fields % elem (k) % mask % nodata) THEN
                  irrigationFlux % mat (i,j) = irrigationFlux % mat (i,j) + &
                                              flux * fields % elem (k) % eta
            END IF
          END DO
        END DO
        
    ELSE
        fields % elem (k) % fluxQ = 0. 
    END IF !doy

END DO

!write diverted and downstream released discharge
IF (time == timeOut) THEN
    WRITE (unitIrrigationDiverted,'(a)', advance = 'no') timeString
    WRITE (unitIrrigationDownstream,'(a)', advance = 'no') timeString
    WRITE (unitIrrigationUpstream,'(a)', advance = 'no') timeString
    
    DO k = 1, fields % count 
   
        r = fields % elem (k) % r
        c = fields % elem (k) % c
        Qdownstream = Qriver % mat (r,c) -  Qirrigation % mat (r,c)
        
        IF ( k < fields % count ) THEN
            WRITE (unitIrrigationDiverted,fmt='(" ",e14.7)', advance = 'no') &
                    fields % elem (k) % fluxQ
            WRITE (unitIrrigationDownstream,fmt='(" ",e14.7)', advance = 'no')&
                    Qdownstream
            WRITE (unitIrrigationUpstream,fmt='(" ",e14.7)', advance = 'no')&
                    Qriver % mat (r,c)
        ELSE
            WRITE (unitIrrigationDiverted,fmt='(" ",e14.7)') &
                    fields % elem (k) % fluxQ
            WRITE (unitIrrigationDownstream,fmt='(" ",e14.7)') Qdownstream
            WRITE (unitIrrigationUpstream,fmt='(" ",e14.7)')  Qriver % mat (r,c)
           
        END IF
    
    END DO
timeOut = timeOut + dtOut

END IF

RETURN

END SUBROUTINE IrrigationUpdate





!==============================================================================
!| Description:
!   populate  array of daily values
FUNCTION SetDailyArray &
!
( value ) &
!
RESULT ( array )

IMPLICIT NONE

!Arguments with intent(in):
CHARACTER (LEN = *) :: value


!Local declarations:
REAL (KIND = float) :: array (365)
REAL (KIND = float) :: scalar
LOGICAL :: error
TYPE (Table) :: valueTable
INTEGER :: i
REAL (KIND = float) :: doyStart, doyStop

!----------------------------end of declarations-------------------------------

!check that value is a number
scalar = StringToFloat (value, error)
IF  ( error ) THEN !value changes in time
    CALL TableNew (value, valueTable)
    array = 0.
    
    DO i = 1, TableGetNrows (valueTable)
       CALL TableGetValue ( valueIn = REAL (i), &
                             tab = valueTable, &
                             keyIn = 'count', &
                             keyOut ='doy-start', &
                             match = 'exact', &
                             valueOut = doyStart )
      
       CALL TableGetValue ( valueIn = REAL (i), &
                             tab = valueTable, &
                             keyIn = 'count', &
                             keyOut ='doy-stop', &
                             match = 'exact', &
                             valueOut = doyStop )
       
       CALL TableGetValue ( valueIn = REAL (i), &
                             tab = valueTable, &
                             keyIn = 'count', &
                             keyOut ='value', &
                             match = 'exact', &
                             valueOut = scalar )
       
       array ( INT (doyStart) : INT (doyStop) ) = scalar
        
    END DO
    
ELSE !no error, value is a scalar
    array = scalar
END IF

RETURN
END FUNCTION SetDailyArray




END MODULE Irrigation