!! Implement models to compute infiltration rate
!|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.1 - 16th June 2016  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 28/May/2014 | Original code |
! | 1.1      | 16/Jun/2016 | added CurveNumber function |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! Implement models to compute infiltration rate:
!
! * [[Philip(function)]] model
!
! * [[GreenAmpt(function)]] model
!
! * [[SCScurveNumber(function)]] model
!
MODULE Infiltration


USE DataTypeSizes, ONLY : &
    ! Imported Type Definitions:
    short, float, double  

USE LogLib,ONLY : &
    ! imported routines:
    Catch 

USE IniLib, ONLY : &
    !Imported types:
    IniList, &
    !Imported routines:
    IniOpen, IniClose , &
    IniReadInt, IniReadReal , &
    KeyIsPresent, IniReadString, &
    IniReadDouble, SectionIsPresent

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

USE StringManipulation, ONLY: &
  !Imported routines:
  ToString

USE GridOperations, ONLY : &
  !imported routines:
  GridByIni, CRSisEqual !, CellArea, &
  !Imported assignment:
  !ASSIGNMENT (=)

USE SoilProperties, ONLY : &
  !Imported definitions:
  SoilType

USE DomainProperties, ONLY : &
    !imported variables:
    mask

USE Richards, ONLY : &
  !imported routines:
  SetRichardsBC,  SolveRichardsBC, &
  !imported variables:
  nsteps, parameters, cell, &
  dx, S, h0grid, drncum, &
  infilcum, evapcum, runoffcum

IMPLICIT NONE 
! Global (i.e. public) Declarations:
INTEGER (KIND = short) :: infiltrationModel

!infiltration model labels
INTEGER (KIND = short), PARAMETER :: SCS_CN = 1 !!Soil Conservation Service Curve Number
INTEGER (KIND = short), PARAMETER :: PHILIPEQ = 2 !!Philip's equation
INTEGER (KIND = short), PARAMETER :: GREEN_AMPT = 3 !!Green Ampt
INTEGER (KIND = short), PARAMETER :: ROSS_BC = 4 !!Richards equations solution by Ross with Broooks and Corey SWRC
INTEGER (KIND = short), PARAMETER :: ROSS_VG = 5 !!Richards equations solution by Ross with Van Genuchten SWRC

!soil
TYPE (SoilType), PRIVATE, ALLOCATABLE :: soils (:)
INTEGER (KIND = short) :: soilDivisions !! number of subdivisions of soil layer

!common properties
TYPE (grid_integer), PRIVATE :: soilTypeMap !read in subroutine SetParametersFromDB
TYPE (grid_real), PUBLIC     :: ksat !!saturated hydraulic conductivity [m/s]
TYPE (grid_real), PUBLIC     :: thetar !!residual water content [m3/m3]
TYPE (grid_real), PUBLIC     :: thetas !!saturated water content [m3/m3]
TYPE (grid_real), PUBLIC     :: psdi !!pore size distribution index by Brooks and Corey [-]
TYPE (grid_real), PUBLIC     :: wiltingPoint !!wilting point [m3/m3]
TYPE (grid_real), PUBLIC     :: fieldCapacity !!field capacity [m3/m3]
TYPE (grid_real), PRIVATE    :: smax !!maximum soil storage  [m] DEBUG non so chi usa questo parametro

!used by PHILIPS, ROSS_BC and ROSS_VG
TYPE (grid_real), PUBLIC     :: psic !!bubbling pressure, air entry value  [m]

!used by SCS_CN
TYPE (grid_real), PUBLIC     :: curveNumber !! curve number [-]
TYPE (grid_real), PUBLIC     :: abstractionRatio !! initial abstraction ratio of the SCS-CN [-] default = 0.2
TYPE (grid_real), PUBLIC     :: storativity !! storativity of the SCS-CN [mm], default = 254
TYPE (grid_real), PUBLIC     :: sEff !! effective soil retention capacity [mm]
TYPE (grid_real), PUBLIC     :: raincum !!cumulated gross rainfall [mm]

!used by GREEN_AMPT
TYPE (grid_real), PUBLIC    :: phy !!wetting front soil suction head  [m]

!used by ROSS_VG
TYPE (grid_real), PRIVATE    :: ptort !!tortuosity index  [-]
TYPE (grid_real), PRIVATE    :: nvg !!n shape parameter for Van Genuchten SWRC  [-]
TYPE (grid_real), PRIVATE    :: mvg !!m shape parameter for Van Genuchten SWRC  [-]
TYPE (grid_real), PRIVATE    :: ksatMatrix !!saturated conductivity of soil "matrix", 
                                !different from ksat when macropores impact. 
                                !Implemented only when vG WRC is used

!used by PHILIPEQ and GREEN_AMPT
TYPE (grid_real), PUBLIC  :: cuminf !!cumulative infiltration [m]
TYPE (grid_real), PUBLIC  :: ism !!soil moisture at the beginning of storm event [m3/m3]



!private declarations:
INTEGER (KIND = short), PRIVATE :: parameterAssigningMethod
!!1 = read each parameter from single file; 2 = assign parameter from soil type map
INTEGER (KIND = short), PRIVATE, PARAMETER :: FROM_FILES = 1
INTEGER (KIND = short), PRIVATE, PARAMETER :: FROM_DB = 2


!public routines:
PUBLIC :: InfiltrationInit
PUBLIC :: Philip
PUBLIC :: GreenAmpt
PUBLIC :: SCScurveNumber 
PUBLIC :: Sat2scn

!private routines:
PRIVATE :: ReadSoilTypes
PRIVATE :: SetParametersFromDB
PRIVATE :: SetParametersFromFile
                
!=======
CONTAINS
!=======
! Define procedures contained in this module. 
  

!==============================================================================
!| Description:
!   Initialize  evapotranspiration computation
SUBROUTINE InfiltrationInit &
!
( inifile, soilMoisture, soilDepth )

IMPLICIT NONE

!Arguments with intent(in):
CHARACTER (LEN = *), INTENT(IN) :: inifile
TYPE (grid_real), INTENT (IN) :: soilMoisture
TYPE (grid_real), INTENT (IN) :: soilDepth

!local declarations:
TYPE (IniList) :: infiltini

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

!open and read configuration file
CALL IniOpen (inifile, infiltini)

!read infiltration model
infiltrationModel = IniReadInt ('model', infiltini)

!read number of subdivisions of soil layer
soilDivisions = IniReadInt ('ross-divisions', infiltini)

!check soil division and infiltration model
IF (soilDivisions > 1 .AND. infiltrationModel == SCS_CN .OR. &
    soilDivisions > 1 .AND. infiltrationModel == PHILIPEQ .OR. &
    soilDivisions > 1 .AND. infiltrationModel == GREEN_AMPT ) THEN
  
     CALL Catch ('warning', 'Infiltration', &
       'Infiltration model requires one subdivision of soil layer, divisions forced to 1')
     
     soilDivisions = 1
END IF

    
!read parameter assigning method 
parameterAssigningMethod = IniReadInt ('parameter-assigning-method', infiltini)

IF (parameterAssigningMethod == FROM_DB) THEN
    !check if all necessary info are defined
    IF (.NOT. KeyIsPresent ('soil-types-file', infiltini) ) THEN
       CALL Catch ('error', 'Infiltration', &
       'missing ''soil-types-file'' keyword in configuration file')
    END IF
    
    IF (.NOT. SectionIsPresent ('soil-type-map', infiltini) ) THEN
       CALL Catch ('error', 'Infiltration', &
       'missing ''soil-type-map'' section in configuration file')
    END IF
    
    !set parameter maps
     CALL SetParametersFromDB (infiltini, infiltrationModel)
    
ELSE
     CALL SetParametersFromFile (infiltini, infiltrationModel)
END IF


!used by Philip model 
IF (infiltrationModel == PHILIPEQ .OR. &
    infiltrationModel == GREEN_AMPT ) THEN
   CALL NewGrid (cuminf, mask, 0.)
   CALL NewGrid (ism, mask, 0.)
END IF

!used by SCS-CN
IF (infiltrationModel == SCS_CN) THEN
   CALL NewGrid (raincum, mask, 0.)
END IF

!when Richards equation solution is required, initialize Ross'parameters
IF (infiltrationModel == ROSS_BC) THEN
  CALL SetRichardsBC (mask, ksat, thetas, thetar, psdi, psic, &
                      soilDepth, soilMoisture, infiltini, soilDivisions)
END IF
    


!  Configuration terminated. Deallocate ini database
CALL IniClose (infiltini)  

RETURN
END SUBROUTINE InfiltrationInit



!==============================================================================
!| Description: 
!  read soil types from external file
SUBROUTINE ReadSoilTypes   & 
  !
  (inifile)           

IMPLICIT NONE

!Arguments with intent(in):
CHARACTER (LEN = *), INTENT(IN) :: inifile !!stores configuration information

!local declarations
TYPE(IniList) :: soilDB
INTEGER (KIND = short) :: nsoils !number of soil types to read
INTEGER (KIND = short) :: i

!------------end of declaration------------------------------------------------ 

!open and read soil types file
CALL IniOpen (inifile, soilDB)

!get number of soil types to read
IF (KeyIsPresent ('soil-types', soilDB) ) THEN !mandatory key
   nsoils = IniReadInt ('soil-types', soilDB)
ELSE
    CALL Catch ('error', 'SoilBalance', &
       'missing ''soil-types'' section in soil types file') 
END IF

!allocate memory and read data
ALLOCATE (soils(nsoils))
DO i = 1, nsoils
  soils (i) % ksat = IniReadDouble ('ksat', soilDB, section = ToString(i)) 
  soils (i) % thetas = IniReadDouble ('thetas', soilDB, section = ToString(i)) 
  soils (i) % thetar = IniReadDouble ('thetar', soilDB, section = ToString(i))
  soils (i) % wp = IniReadDouble ('wp', soilDB, section = ToString(i))
  soils (i) % fc = IniReadDouble ('fc', soilDB, section = ToString(i))
  soils (i) % psic = IniReadDouble ('psic', soilDB, section = ToString(i))
  soils (i) % psdi = IniReadDouble ('psdi', soilDB, section = ToString(i))
  soils (i) % phy = IniReadDouble ('phy', soilDB, section = ToString(i))
  soils (i) % smax = IniReadDouble ('smax', soilDB, section = ToString(i))
  soils (i) % m = IniReadDouble ('m', soilDB, section = ToString(i))
  soils (i) % n = IniReadDouble ('n', soilDB, section = ToString(i))
  soils (i) % pp = IniReadDouble ('p', soilDB, section = ToString(i))
  soils (i) % kx = IniReadDouble ('ksat-matrix', soilDB, section = ToString(i))
  soils (i) % c = IniReadDouble ('c', soilDB, section = ToString(i))
  soils (i) % s0 = IniReadDouble ('s0', soilDB, section = ToString(i))
  soils (i) % cn = IniReadDouble ('cn', soilDB, section = ToString(i))
  
END DO

!deallocate soilDB
CALL IniClose (soilDB)

RETURN
END SUBROUTINE ReadSoilTypes
  

!==============================================================================
!| Description: 
!   set parameter maps from soil database
SUBROUTINE SetParametersFromDB   & 
!
  (iniDB, model)

IMPLICIT NONE

! Arguments with intent(in):

TYPE (IniList), INTENT(IN) :: iniDB
INTEGER (KIND = short), INTENT(IN) :: model !infiltration model

!Local declaration:
CHARACTER (LEN = 1000) :: soilTypeFile
INTEGER (KIND = short) ::i, j, id

!------------end of declaration------------------------------------------------ 
CALL Catch ('info', 'SoilBalance: ', 'setting soil parameters from database: ', &
          argument = IniReadString('soil-types-file', iniDB) )

!load soil types
soilTypeFile = IniReadString('soil-types-file', iniDB)
CALL ReadSoilTypes (soilTypeFile)


! load soil type map
CALL GridByIni (iniDB, soilTypeMap, section = 'soil-type-map')


!Set parameter maps used by all models

!first allocate memory
CALL NewGrid (ksat, soilTypeMap) 
CALL NewGrid (thetar, soilTypeMap) 
CALL NewGrid (thetas, soilTypeMap)
CALL NewGrid (wiltingPoint, soilTypeMap)
CALL NewGrid (fieldCapacity, soilTypeMap)
CALL NewGrid (psdi, soilTypeMap)

!then assigh parameters
DO i = 1, soilTypeMap % idim
  DO j = 1, soilTypeMap % jdim
     IF (soilTypeMap % mat (i,j) /= soilTypeMap % nodata) THEN
       id = soilTypeMap % mat (i,j)
       ksat % mat (i,j) = soils (id) % ksat 
       thetar % mat (i,j) = soils (id) % thetar
       thetas % mat (i,j) = soils (id) % thetas
       wiltingPoint % mat (i,j) = soils (id) % wp
       fieldCapacity % mat (i,j) = soils (id) % fc
       psdi % mat (i,j) = soils (id) % psdi
     END IF
  END DO
END DO  

IF (model == SCS_CN) THEN !read supplementary parameters required by Curve Number 
  !first allocate memory
  CALL NewGrid (curveNumber, soilTypeMap) 
  CALL NewGrid (abstractionRatio, soilTypeMap)
  CALL NewGrid (storativity, soilTypeMap)
  
  !then assign parameters
  DO i = 1, soilTypeMap % idim
    DO j = 1, soilTypeMap % jdim
       IF (soilTypeMap % mat (i,j) /= soilTypeMap % nodata) THEN
         id = soilTypeMap % mat (i,j)
         curveNumber % mat (i,j) = soils (id) % cn 
         abstractionRatio % mat (i,j) = soils (id) % c
         storativity % mat (i,j) = soils (id) % s0
       END IF
    END DO
  END DO
END IF

IF (model == PHILIPEQ) THEN !read supplementary parameters required by Philips 
  !first allocate memory
  CALL NewGrid (psic, soilTypeMap) 
  
  !then assigh parameters
  DO i = 1, soilTypeMap % idim
    DO j = 1, soilTypeMap % jdim
       IF (soilTypeMap % mat (i,j) /= soilTypeMap % nodata) THEN
         id = soilTypeMap % mat (i,j)
         psic % mat (i,j) = soils (id) % psic 
       END IF
    END DO
  END DO
END IF

IF (model == GREEN_AMPT) THEN !read supplementary parameters required by Green Ampt 
  !first allocate memory
  CALL NewGrid (phy, soilTypeMap) 
  
  !then assigh parameters
  DO i = 1, soilTypeMap % idim
    DO j = 1, soilTypeMap % jdim
       IF (soilTypeMap % mat (i,j) /= soilTypeMap % nodata) THEN
         id = soilTypeMap % mat (i,j)
         phy % mat (i,j) = soils (id) % phy 
       END IF
    END DO
  END DO
END IF

IF (model == ROSS_BC) THEN !read supplementary parameters required by Ross Brooks and Corey 
  !first allocate memory
  CALL NewGrid (psic, soilTypeMap) 
  
  !then assigh parameters
  DO i = 1, soilTypeMap % idim
    DO j = 1, soilTypeMap % jdim
       IF (soilTypeMap % mat (i,j) /= soilTypeMap % nodata) THEN
         id = soilTypeMap % mat (i,j)
         psic % mat (i,j) = soils (id) % psic 
       END IF
    END DO
  END DO
END IF

IF (model == ROSS_VG) THEN !read supplementary parameters required by Ross Van Genuchten 
  !first allocate memory
  CALL NewGrid (psic, soilTypeMap) 
  CALL NewGrid (nvg, soilTypeMap) 
  CALL NewGrid (mvg, soilTypeMap) 
  CALL NewGrid (ptort, soilTypeMap)
  CALL NewGrid (ksatMatrix, soilTypeMap) 
  
  !then assigh parameters
  DO i = 1, soilTypeMap % idim
    DO j = 1, soilTypeMap % jdim
       IF (soilTypeMap % mat (i,j) /= soilTypeMap % nodata) THEN
         id = soilTypeMap % mat (i,j)
         psic % mat (i,j) = soils (id) % psic 
         nvg % mat (i,j) = soils (id) % n 
         mvg % mat (i,j) = soils (id) % m 
         mvg % mat (i,j) = soils (id) % m 
         ksatMatrix % mat (i,j) = soils (id) % kx
       END IF
    END DO
  END DO
END IF

RETURN
END SUBROUTINE SetParametersFromDB

  
  
  
!==============================================================================
!| Description: 
!   set parameter maps from separate files
SUBROUTINE SetParametersFromFile   & 
!
  (iniDB, model)

IMPLICIT NONE

! Arguments with intent(in):

TYPE (IniList), INTENT(IN) :: iniDB
INTEGER (KIND = short), INTENT(IN) :: model !infiltration model

!Local declaration:

!------------end of declaration------------------------------------------------ 

!Set parameter maps used by all models

!Saturated Conductivity
IF (SectionIsPresent('conductivity', iniDB)) THEN
    CALL GridByIni (iniDB, ksat, section = 'conductivity')
    CALL Catch ('info', 'SoilBalance: ', 'loading saturated conductivity from file: ', &
          argument = IniReadString('file', iniDB, section = 'conductivity'))
    IF  ( .NOT. CRSisEqual (mask = mask, grid = ksat, checkCells = .TRUE.) ) THEN
        CALL Catch ('error', 'SoilBalance',   &
			    'wrong spatial reference in conductivity map' )
    END IF
ELSE
    CALL Catch ('error', 'SoilBalance: ',   &
			    'missing conductivity section in configuration file' )
END IF
 
!Saturated water content
IF (SectionIsPresent('saturated-water-content', iniDB)) THEN
    CALL GridByIni (iniDB, thetas, section = 'saturated-water-content')
    CALL Catch ('info', 'SoilBalance: ', 'loading saturated water content from file: ', &
          argument = IniReadString('file', iniDB, section = 'saturated-water-content'))
    IF  ( .NOT. CRSisEqual (mask = mask, grid = thetas, checkCells = .TRUE.) ) THEN
        CALL Catch ('error', 'SoilBalance',   &
			    'wrong spatial reference in saturated-water-content' )
    END IF
ELSE
    CALL Catch ('error', 'SoilBalance: ',   &
			    'missing saturated-water-content section in configuration file' )
END IF

!Residual water content
IF (SectionIsPresent('residual-water-content', iniDB)) THEN
    CALL GridByIni (iniDB, thetar, section = 'residual-water-content')
    CALL Catch ('info', 'SoilrBalance: ', 'loading residual water content from file: ', &
          argument = IniReadString('file', iniDB, section = 'residual-water-content'))
    IF  ( .NOT. CRSisEqual (mask = mask, grid = thetar, checkCells = .TRUE.) ) THEN
        CALL Catch ('error', 'SoilBalance',   &
			    'wrong spatial reference in residual-water-content' )
    END IF
ELSE
    CALL Catch ('error', 'SoilWaterBalance: ',   &
			    'missing residual-water-content section in configuration file' )
END IF


!Wilting Point
IF (SectionIsPresent('wilting-point', iniDB)) THEN
    CALL GridByIni (iniDB, wiltingPoint, section = 'wilting-point')
    CALL Catch ('info', 'SoilBalance: ', 'loading wilting point from file: ', &
          argument = IniReadString('file', iniDB, section = 'wilting-point'))
    IF  ( .NOT. CRSisEqual (mask = mask, grid = wiltingPoint, checkCells = .TRUE.) ) THEN
        CALL Catch ('error', 'SoilBalance',   &
			    'wrong spatial reference in wilting-point' )
    END IF
ELSE
    CALL Catch ('error', 'SoilBalance: ',   &
			    'missing wilting-point section in configuration file' )
END IF
    
!Field Capacity
IF (SectionIsPresent('field-capacity', iniDB)) THEN
    CALL GridByIni (iniDB, fieldCapacity, section = 'field-capacity')
    CALL Catch ('info', 'SoilBalance: ', 'loading field capacity from file: ', &
          argument = IniReadString('file', iniDB, section = 'field-capacity'))
    IF  ( .NOT. CRSisEqual (mask = mask, grid = fieldCapacity, checkCells = .TRUE.) ) THEN
        CALL Catch ('error', 'SoilBalance',   &
			    'wrong spatial reference in field-capacity' )
    END IF
ELSE
    CALL Catch ('error', 'SoilBalance: ',   &
			    'missing field-capacity section in configuration file' )
END IF

!Pore size distribution index
IF (SectionIsPresent('pore-size-index', iniDB)) THEN
    CALL GridByIni (iniDB, psdi, section = 'pore-size-index')
    CALL Catch ('info', 'SoilBalance: ', 'loading pore size index from file: ', &
          argument = IniReadString('file', iniDB, section = 'pore-size-index'))
    IF  ( .NOT. CRSisEqual (mask = mask, grid = psdi, checkCells = .TRUE.) ) THEN
        CALL Catch ('error', 'SoilBalance',   &
			    'wrong spatial reference in pore-size-index')
    END IF
ELSE
    CALL Catch ('error', 'SoilBalance: ',   &
			    'missing pore-size-idex section in configuration file' )
END IF

!bubbling pressure, air entry value
IF (SectionIsPresent('bubble-pressure', iniDB)) THEN
    CALL GridByIni (iniDB, psic, section = 'bubble-pressure')
    CALL Catch ('info', 'SoilBalance: ', 'loading bubbling pressure from file: ', &
            argument = IniReadString('file', iniDB, section = 'bubble-pressure'))
    IF  ( .NOT. CRSisEqual (mask = mask, grid = psic, checkCells = .TRUE.) ) THEN
        CALL Catch ('error', 'SoilBalance',   &
			    'wrong spatial reference in bubbling pressure')
    END IF
ELSE
      CALL Catch ('error', 'SoilBalance: ',   &
			      'missing bubble-pressure section in configuration file' )
END IF


IF (model == SCS_CN) THEN !read supplementary parameters required by Curve Number 
  
  !Curve Number
	IF (SectionIsPresent('curve-number', iniDB)) THEN
        CALL GridByIni (iniDB, curveNumber, section = 'curve-number')
        CALL Catch ('info', 'SoilBalance: ', 'loading curve number from file: ', &
             argument = IniReadString('file', iniDB, section = 'curve-number'))
        IF  ( .NOT. CRSisEqual (mask = mask, grid = curveNumber, checkCells = .TRUE.) ) THEN
           CALL Catch ('error', 'SoilBalance',   &
			        'wrong spatial reference in curve number')
        END IF
  ELSE
      CALL Catch ('error', 'SoilBalance: ',   &
			      'missing curve-number section in configuration file' )
  END IF
  
  
  !Initial abstraction ratio
	IF (SectionIsPresent('abstraction-ratio', iniDB)) THEN
        CALL GridByIni (iniDB, abstractionRatio, section = 'abstraction-ratio')
        CALL Catch ('info', 'SoilBalance: ', 'loading initial abstraction ratio from file: ', &
             argument = IniReadString('file', iniDB, section = 'abstraction-ratio'))
        IF  ( .NOT. CRSisEqual (mask = mask, grid = abstractionRatio, checkCells = .TRUE.) ) THEN
           CALL Catch ('error', 'SoilBalance',   &
			        'wrong spatial reference in initial abstraction ratio')
        END IF
  ELSE
      CALL Catch ('error', 'SoilBalance: ',   &
			      'missing abstraction-ratio section in configuration file' )
  END IF
  
  !SCS-CN storativity
	IF (SectionIsPresent('storativity', iniDB)) THEN
        CALL GridByIni (iniDB, storativity, section = 'storativity')
        CALL Catch ('info', 'SoilBalance: ', 'loading storativity from file: ', &
             argument = IniReadString('file', iniDB, section = 'storativity'))
        IF  ( .NOT. CRSisEqual (mask = mask, grid = storativity, checkCells = .TRUE.) ) THEN
           CALL Catch ('error', 'SoilBalance',   &
			        'wrong spatial reference in storativity')
        END IF
  ELSE
      CALL Catch ('error', 'SoilBalance: ',   &
			      'missing storativity section in configuration file' )
  END IF
  
END IF  

!IF (model == PHILIPEQ) THEN !read supplementary parameters required by Philips 
!  
!  !bubbling pressure, air entry value
!	IF (SectionIsPresent('bubble-pressure', iniDB)) THEN
!        CALL GridByIni (iniDB, psic, section = 'bubble-pressure')
!        CALL Catch ('info', 'SoilBalance: ', 'loading bubbling pressure from file: ', &
!             argument = IniReadString('file', iniDB, section = 'bubble-pressure'))
!        IF  ( .NOT. CRSisEqual (mask = mask, grid = psic, checkCells = .TRUE.) ) THEN
!           CALL Catch ('error', 'SoilBalance',   &
!			        'wrong spatial reference in bubbling pressure')
!        END IF
!  ELSE
!      CALL Catch ('error', 'SoilBalance: ',   &
!			      'missing bubble-pressure section in configuration file' )
!  END IF
!
!END IF

IF (model == GREEN_AMPT) THEN !read supplementary parameters required by Green Ampt 
  
  !wetting front soil suction head
	IF (SectionIsPresent('front-suction-head', iniDB)) THEN
        CALL GridByIni (iniDB, phy, section = 'front-suction-head')
        CALL Catch ('info', 'SoilBalance: ', 'loading wetting front suction head from file: ', &
             argument = IniReadString('file', iniDB, section = 'front-suction-head'))
        IF  ( .NOT. CRSisEqual (mask = mask, grid = phy, checkCells = .TRUE.) ) THEN
           CALL Catch ('error', 'SoilBalance',   &
			        'wrong spatial reference in wetting front suction head')
        END IF
  ELSE
      CALL Catch ('error', 'SoilBalance: ',   &
			      'missing front-suction-head section in configuration file' )
  END IF
  
END IF

IF (model == ROSS_BC) THEN !read supplementary parameters required by Ross Brooks and Corey 
  
  !bubbling pressure, air entry value
	IF (SectionIsPresent('bubble-pressure', iniDB)) THEN
        CALL GridByIni (iniDB, psic, section = 'bubble-pressure')
        CALL Catch ('info', 'SoilBalance: ', 'loading bubbling pressure from file: ', &
             argument = IniReadString('file', iniDB, section = 'bubble-pressure'))
        IF  ( .NOT. CRSisEqual (mask = mask, grid = psic, checkCells = .TRUE.) ) THEN
           CALL Catch ('error', 'SoilBalance',   &
			        'wrong spatial reference in bubbling pressure')
        END IF
  ELSE
      CALL Catch ('error', 'SoilBalance: ',   &
			      'missing bubble-pressure section in configuration file' )
  END IF

  
END IF

IF (model == ROSS_VG) THEN !read supplementary parameters required by Ross Van Genuchten 
  
  !bubbling pressure, air entry value
	IF (SectionIsPresent('bubble-pressure', iniDB)) THEN
        CALL GridByIni (iniDB, psic, section = 'bubble-pressure')
        CALL Catch ('info', 'SoilBalance: ', 'loading bubbling pressure from file: ', &
             argument = IniReadString('file', iniDB, section = 'bubble-pressure'))
        IF  ( .NOT. CRSisEqual (mask = mask, grid = psic, checkCells = .TRUE.) ) THEN
           CALL Catch ('error', 'SoilBalance',   &
			        'wrong spatial reference in bubbling pressure')
        END IF
  ELSE
      CALL Catch ('error', 'SoilBalance: ',   &
			      'missing bubble-pressure section in configuration file' )
  END IF
  
  !n shape parameter for Van Genuchten SWRC
	IF (SectionIsPresent('n-van-genuchten', iniDB)) THEN
        CALL GridByIni (iniDB, nvg, section = 'n-van-genuchten')
        CALL Catch ('info', 'SoilBalance: ', 'loading n shape parameter from file: ', &
             argument = IniReadString('file', iniDB, section = 'n-van-genuchten'))
        IF  ( .NOT. CRSisEqual (mask = mask, grid = nvg, checkCells = .TRUE.) ) THEN
           CALL Catch ('error', 'SoilBalance',   &
			        'wrong spatial reference in n shape parameter')
        END IF
  ELSE
      CALL Catch ('error', 'SoilBalance: ',   &
			      'missing n-van-genuchten section in configuration file' )
  END IF
  
  !m shape parameter for Van Genuchten SWRC
	IF (SectionIsPresent('m-van-genuchten', iniDB)) THEN
        CALL GridByIni (iniDB, mvg, section = 'm-van-genuchten')
        CALL Catch ('info', 'SoilBalance: ', 'loading m shape parameter from file: ', &
             argument = IniReadString('file', iniDB, section = 'm-van-genuchten'))
        IF  ( .NOT. CRSisEqual (mask = mask, grid = mvg, checkCells = .TRUE.) ) THEN
           CALL Catch ('error', 'SoilBalance',   &
			        'wrong spatial reference in m shape parameter')
        END IF
  ELSE
      CALL Catch ('error', 'SoilBalance: ',   &
			      'missing m-van-genuchten section in configuration file' )
  END IF
  
  !tortuosity index
	IF (SectionIsPresent('tortuosity-index', iniDB)) THEN
        CALL GridByIni (iniDB, ptort, section = 'tortuosity-index')
        CALL Catch ('info', 'SoilBalance: ', 'loading tortuosity index from file: ', &
             argument = IniReadString('file', iniDB, section = 'tortuosity-index'))
        IF  ( .NOT. CRSisEqual (mask = mask, grid = ptort, checkCells = .TRUE.) ) THEN
           CALL Catch ('error', 'SoilBalance',   &
			        'wrong spatial reference in tortuosity index')
        END IF
  ELSE
      CALL Catch ('error', 'SoilBalance: ',   &
			      'missing tortuosity-index section in configuration file' )
  END IF
  
  !saturated conductivity of soil "matrix"
	IF (SectionIsPresent('conductivity-matrix', iniDB)) THEN
        CALL GridByIni (iniDB, ksatMatrix, section = 'conductivity-matrix')
        CALL Catch ('info', 'SoilBalance: ', 'loading soil matrix conductivity from file: ', &
             argument = IniReadString('file', iniDB, section = 'conductivity-matrix'))
        IF  ( .NOT. CRSisEqual (mask = mask, grid = ksatMatrix, checkCells = .TRUE.) ) THEN
           CALL Catch ('error', 'SoilBalance',   &
			        'wrong spatial reference in soil matrix conductivity')
        END IF
  ELSE
      CALL Catch ('error', 'SoilBalance: ',   &
			      'missing conductivity-matrix section in configuration file' )
  END IF

END IF

RETURN
  END SUBROUTINE SetParametersFromFile  





  
  
!==============================================================================
!| calculates the actual rate of infiltration (m/s) according to 
! Curve Number model modified for continuous simulation
!
! References:
!
!  Ravazzani, G., Mancini, M., Giudici, I., & Amadio, P.. (2007). 
!     Effects of soil moisture parameterization on a real- time flood 
!     forecasting system based on rainfall thresholds. 
!     IAHS publication, 313, 407–416.
!
!  Rabuffetti, D., Ravazzani, G., Corbari, C., & Mancini, M.. (2008). 
!     Verification of operational quantitative discharge forecast (QDF) 
!     for a regional warning system – the AMPHORE case studies in the 
!     upper Po river. Natural hazards and earth system sciences, 8, 161–173.
!
!  Ravazzani, G., Amengual, A., Ceppi, A., Homar, V., Romero, R.,
!     Lombardi, G., & Mancini, M.. (2016). Potentialities of ensemble 
!     strategies for flood forecasting over the Milano urban area. 
!     Journal of hydrology, 539, 237-253.
!
FUNCTION SCScurveNumber &
!
(rain, raincum, c, s, dt, runoff) &
!
RESULT (inf)


IMPLICIT NONE

!Arguments with intent in
REAL (KIND = float), INTENT(IN) :: rain !!rainfall rate (m/s)
REAL (KIND = float), INTENT(IN) :: c !!initial abstraction ratio (-)
REAL (KIND = float), INTENT(IN) :: s !!soil retention capacity (mm)
INTEGER (KIND = short), INTENT(IN) :: dt !!time step (s)

!Arguments with intent inout:
REAL (KIND = float), INTENT(INOUT) :: raincum !!cumulated rainfall from start of storm (mm)
                                              !! is used and updated to be ready for the 
                                              !! following step

!Arguments with intent out
REAL (KIND = double),  INTENT(OUT) :: runoff !!runoff rate (m/s)

!local declarations:
REAL (KIND = float) :: inf !!infiltration rate (m/s)
REAL (KIND = float) :: raincump  !!cumulative precipitation of previous time step (mm)
REAL (KIND = float) :: runoffp  !!cumulative runoff of previous time step (mm)


!-------------------------end of declarations----------------------------------
!if retention capacity = 0 (soil is saturated)
! rainfall is transformed to runoff 
IF (s == 0.) THEN
		raincum = raincum + rain * 1000. * dt !update cumulative rainfall
		runoff = rain
		inf = 0.
		RETURN !no need to continue
END IF

!calculate runoff and infiltration rate  

!update cumulative precipitation
raincump = raincum !save cumulated precipitation amount of previous time step
raincum = raincum + rain * dt * 1000. !(mm)
  
!runoff at current time
IF(raincum >= c*s) THEN
	runoff = ((raincum - c*s )**2.) / (raincum + (1.-c) * s)
ELSE
	runoff = 0.
END IF
     
!runoff at previous time
IF(raincump >= c*S) THEN
	runoffp = ((raincump - c*S )**2.) / (raincump + (1.-c) * s)
ELSE
	runoffp = 0.
END IF
    
!actual runoff in m/s
runoff = (runoff - runoffp) / 1000. / dt
    
!actual infiltration rate in m/s
inf = rain - runoff

RETURN
END FUNCTION SCScurveNumber

!==============================================================================
!| compute actual soil retention capacity S of the SCS-CN method  
! given soil saturation
SUBROUTINE Sat2scn &
  !
  (cn2, s0, soilsat, scn)

IMPLICIT NONE

!arguments with intent in:
REAL (KIND = float), INTENT(IN) :: cn2 !!curve number for mean soil moisture 
REAL (KIND = float), INTENT(IN) :: s0 !!storativity (mm)
REAL (KIND = float), INTENT(IN) :: soilsat !!soil saturation

!arguments with intent out:
REAL (KIND = float), INTENT(OUT) :: scn !!soil retention capacity [mm]

!local declarations:
REAL (KIND = float) :: cn1, cn3, s1, s2, s3
!-------------------------end of declarations----------------------------------
	!compute cn1 and cn3 from cn2
	cn1 = cn2 - (20. * (100. - cn2) ) / &
        (100. - cn2 + EXP(2.533 - .0636 * (100. - cn2) ) )
  
	IF (cn1 < 0.4 * CN2) cn1 = .4 * cn2
  
	cn3 = cn2 * EXP( .00673 * (100. - cn2) )
  
	!compute s1, s2, and s3
	s1 = s0 * (100. / cn1 - 1.)
	s2 = s0 * (100. / cn2 - 1.)
	s3 = s0 * (100. / cn3 - 1.)
	
	!compute actual S
	scn = s1 - soilsat * (s1 - s3)    

RETURN
  END SUBROUTINE sat2scn

!==============================================================================
!| calculates the actual rate of infiltration (m/s) according to 
! Philip's equation. Use time compression approximation.
!
! References:
!
! Philip, J. R. 1954. An infiltration equation with physical 
!  significance. Soil Science, 77: 153-157. 
!
! Milly, P. C. D., An event-based simulation model of moisture and 
!  energy fluxes at a  bare soil surface, 
!  Water Resour. Res., 22, 1680-! 692, 1986
!
! Sivapalan, M., Beven, K., and Wood, E. F. (1987), On hydrologic similarity: 
!   2 . A scaled model of storm runoff production, Water Resour. 
!   Res., 23( 12), 2266– 2278, doi:10.1029/WR023i012p02266.
!
FUNCTION Philip &
!
(rain, psic, ksat, theta, thetas, thetar, psdi, itheta, dt, cuminf) &
!
RESULT (inf)


IMPLICIT NONE

!Arguments with intent in
REAL (KIND = float), INTENT(IN) :: rain !!rainfall rate (m/s)
REAL (KIND = float), INTENT(IN) :: psic !!bubbling pressure (m)
REAL (KIND = float), INTENT(IN) :: ksat !!saturated hydraulic conductivity (m/s)
REAL (KIND = float), INTENT(IN) :: theta !!volumetric water content (m3/m3)
REAL (KIND = float), INTENT(IN) :: thetas !!saturated volumetric water content (m3/m3)
REAL (KIND = float), INTENT(IN) :: thetar !!residual volumetric water content (m3/m3)
REAL (KIND = float), INTENT(IN) :: psdi !!Brooks & Corey pore size distribution index (-)
INTEGER (KIND = short), INTENT(IN) :: dt !!time step (s)
REAL (KIND = float), INTENT(IN)  :: itheta !!initial soil moisture

!Arguments with intent inout:
REAL (KIND = float), INTENT(INOUT) :: cuminf !!cumulative infiltration (m)

!local declarations:
REAL (KIND = float)  :: inf !!infiltration rate (m/s)
REAL (KIND = float)  :: sorp  !!sorptivity
REAL (KIND = float)  :: cc = 0. !!gravity term in Philip's infiltration equation
REAL (KIND = float)  :: deltrz  !! difference between soil moisture at beginning
                                !! of storm event and residual soil moisture
REAL (KIND = float) :: bcgamm !!2.0 + 3.0 * psdi
REAL (KIND = float) :: infcap !!infiltration capacity (m/s)


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

  ! Sivapalan et al. (1987)
  bcgamm = 2.0 + 3.0 * psdi
  sorp = (((2. * ksat *( (thetas-itheta)**2.) * psic ) / &
         (thetas-thetar))*((1./(bcgamm + 0.5 * psdi - 1.)) + &
         ((thetas-thetar)/ (thetas-itheta) ) ) ) ** 0.5

  deltrz = itheta - thetar
  IF (deltrz <= 0.) deltrz = 0.
  cc = 0.5 * (1. + ( (deltrz / (thetas - thetar))**(bcgamm / psdi) ) )


!If soil is saturated then set infiltration to zero
IF (theta >= thetas) THEN
  inf = 0.
  
!If surface is unsaturated then calculate infiltration  
ELSE
  !If first step with infiltration then all precipitation is
  !infiltrated to avoid division by zero
  IF (cuminf == 0.) THEN
    inf = rain 

    
  !calculate infiltration capacity as a function of cumulative
  ! infiltration for all other time steps of storm (Milly, 1986)  
  ELSE
     infcap = cc * ksat * (1. + (1./ ((( 1. + (( 4. * cc * ksat * cuminf ) / &
              (sorp**2.)))**0.5) - 1.))) 
     !Take the actual infiltration rate as the minimum of the
     ! precipitation rate or the infiltration capacity.
     inf = MIN (infcap,rain)
  END IF
  
  
  !update cumulative infiltration
  cuminf = cuminf + inf * dt

END IF

RETURN

END FUNCTION Philip



!==============================================================================
!| calculates the actual rate of infiltration (m/s) according to 
! Green-Ampt equation for unsteady rainfall.
!
! References:
!
!   Green, WH., and Ampt, GA., 1911. Studies on soil physics I. 
!      Flow of air and water through soils. J. Agric. Sci. 4: 1- 24.
!
!   Chow, V., Maidment, D., Mays, L. (1988) Applied Hydrology,
!      Surface water, pp. 142-145, McGraw-Hill
!
FUNCTION GreenAmpt &
!
(storm, rain,  ksat, theta, thetas, thetar, phy, itheta, dt, cuminf) &
!
RESULT (inf)


IMPLICIT NONE

!Arguments with intent in
INTEGER (KIND = short), INTENT(IN) :: storm !!when = 1 first steo of storm event
REAL (KIND = float), INTENT(IN) :: rain !!rainfall rate (m/s)
REAL (KIND = float), INTENT(IN) :: ksat !!saturated hydraulic conductivity (m/s)
REAL (KIND = float), INTENT(IN) :: theta !!volumetric water content (m3/m3)
REAL (KIND = float), INTENT(IN) :: thetas !!saturated volumetric water content (m3/m3)
REAL (KIND = float), INTENT(IN) :: thetar !!residual volumetric water content (m3/m3)
REAL (KIND = float), INTENT(IN) :: phy !!suction head across the wetting front (m)
INTEGER (KIND = short), INTENT(IN) :: dt !!time step (s)
REAL (KIND = float), INTENT(IN)  :: itheta !!initial soil moisture

!Arguments with intent inout:
REAL (KIND = float), INTENT(INOUT) :: cuminf !!cumulative infiltration (m)


!local declarations:
REAL (KIND = float) :: inf !!infiltration rate (m/s)
REAL (KIND = float)  :: deltrz  !! difference between soil moisture at beginning
                                !! of storm event and residual soil moisture
REAL (KIND = float) :: test  !!test value for the cumulative infiltration (m)
REAL (KIND = float) :: ftest !!test value for the cumulative infiltration (m)
INTEGER (KIND = short) :: i !loop
INTEGER (KIND = short) :: maxiter = 100 !maximum number of iteration for 
                                        ! computing cumulative infiltration
REAL (KIND = float) :: tol = 0.00001 !tolerance to compute cumulative infiltration (m)
REAL (KIND = float) :: deltatfirst !!for ponding time (case 3)
REAL (KIND = float) :: Fponding !!cumulative infiltration at ponding (case 3)
REAL (KIND = float) :: deltatpond !! new deltat calculated for ponding (case 3)



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

!if this is the first step of storm event:  
! initialize variables
IF ( storm == 1) THEN

  deltrz = thetas - theta
  IF (deltrz <= 0.) THEN 
    deltrz = 0.
  END IF 
  
  !assume all precipitation is infiltrated
  cuminf = rain * dt
  inf = rain
  RETURN
END IF

!If soil is saturated then set infiltration to zero
IF (theta >= thetas) THEN
  inf = 0.
  RETURN
!If surface is unsaturated then calculate infiltration  
ELSE
  
  IF (storm == 2) THEN !storm is running 
      
     !calculate potential infiltration rate at the beginning of the interval
     !from the known value of cumulative infiltration (eq 5.4.1)
      inf = ksat * (phy * (thetas-itheta) / cuminf + 1.)
      
      !compare potential infiltration to rainfall rate
      IF (inf <= rain) THEN !ponding occurs throughout the interval (case 1)
      
         test = Ksat*dt !first guess for cumulative infiltration
         DO i = 1, maxiter
             IF (i == maxiter) THEN
                CALL Catch ('warning', 'Infiltration', 'maximum number of iterations &
                             reached while updating cumulative infiltration &
                             with green-ampt case 1' )
             END IF
             ftest = cuminf + ksat*dt + (phy*(thetas-itheta))* &
                     LOG ((test+(phy*(thetas-itheta)))/(cuminf+(phy*(thetas-itheta))))
 
             IF (ABS(ftest-test) <= tol) THEN 
                cuminf = ftest
                !compute infiltration rate from cumulative infiltration
                inf = ksat * (phy * (thetas-itheta) / cuminf + 1.)
                IF (inf > rain) inf = rain
                EXIT
             ELSE 
               test = ftest
             END IF 
         END DO 
               
      ELSE !no ponding at the beginning of the interval
      
        !calculate tentative values Assuming that all the rain
        !water is penetrated into the soil      
        test = cuminf + rain * dt
        inf = ksat * (phy * (thetas-itheta) / test + 1.)
        
        IF (inf <= rain) THEN !ponding occurs during the interval (case 3)
           Fponding = ksat * phy * (thetas-itheta) / (rain - ksat)
           deltatfirst = (Fponding - cuminf) / rain
           deltatpond = dt - deltatfirst
           
           !compute cumulative infiltration by substituting cuminf = Fponding
           !and dt = deltatpond
           test = Fponding !Ksat*dt !first guess for cumulative infiltration
           DO i = 1, maxiter
               IF (i == maxiter) THEN
                 CALL Catch ('warning', 'Infiltration', 'maximum number of iterations &
                             reached while updating cumulative infiltration &
                              with green-ampt case 3' )
               END IF
                 ftest = Fponding + ksat*deltatpond + (phy*(thetas-itheta))* &
                     LOG ((test+(phy*(thetas-itheta)))/(Fponding+(phy*(thetas-itheta))))
 
               IF (ABS(ftest-test) <= tol) THEN 
                  cuminf = ftest
                  !compute infiltration rate from cumulative infiltration
                  inf = ksat * (phy * (thetas-itheta) / cuminf + 1.)
                  IF (inf > rain) inf = rain
                  EXIT
               ELSE 
                 test = ftest
               END IF 
            END DO
            
               
        ELSE !no ponding throghout the interval (case 2)
           cuminf = test
           inf = rain
        END IF
        
      END IF
  
  ELSE
     inf = 0.
  END IF
END IF  

END FUNCTION GreenAmpt 



END MODULE Infiltration