!! Simulate quasi 3D groundwater flux and river-groundwater interaction
!|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 - 12th February 2023 
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 14/Sep/2011 | Original code |
! | 1.0      | 22/Feb/2023 | Code rewritten |
  
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
!Simulate quasi 3D groundwater flux and river-groundwater interaction
! 
!### References
!
! Ravazzani, G., Rametta, D., & Mancini, M.. (2011). Macroscopic cellular 
!   automata for groundwater modelling: a first approach. Environmental 
!   modelling & software, 26(5), 634–643.
!
! Ravazzani, G., Curti, D., Gattinoni, P., Della Valentina, S., 
!   Fiorucci, A., & Rosso, R.. (2016). Assessing groundwater 
!   contribution to streamflow of a large alpine river with heat 
!   tracer methods and hydrological modelling. River research 
!   and applications, 32(5), 871-884.
!
! Ravazzani, G., Barbero, S., Salandin, A., Senatore, A.,
!  & Mancini, M.. (2015). An integrated hydrological model 
!  for assessing climate change impacts on water resources 
!  of the upper po river basin. Water resources management, 
!  29 (4), 1193-1215.
!
MODULE Groundwater

! Modules used:

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

USE DomainProperties, ONLY : &
    !imported variables:
    mask

USE GridLib, ONLY : &
    !imported Type definitions:
    grid_real, grid_integer, &
    !imported routines:
    NewGrid, GridDestroy, &
    !Imported parameters:
    NET_CDF

USE GridOperations, ONLY : &
  !Imported routines:
  GridByIni, CellArea,  &
  !Imported operators:
  ASSIGNMENT (=)

USE IniLib, ONLY : &
  !Imported definitions:
  IniList, &
  !Imported routines:
  IniOpen, IniClose, &
  SectionIsPresent, KeyIsPresent, &
  SubSectionIsPresent, &
  IniReadReal, IniReadInt

USE LogLib, ONLY: &
  !Imported routines:
  Catch

USE StringManipulation, ONLY : &
  !Imported routines
  ToString

USE Chronos, ONLY : &
  !Imported definitions:
  DateTime , &
  !Imported variables:
  timeString, &
 !Imported operations:
  OPERATOR (+), &
  OPERATOR (==), &
  ASSIGNMENT (=)

USE ObservationalNetworks, ONLY : &
  !Imported routines:
  ReadMetadata, CopyNetwork, &
  WriteMetadata, DestroyNetwork, &
  AssignDataFromGrid, WriteData, &
  !Imported defined variable:
  ObservationalNetwork

USE Utilities, ONLY : &
  !Imported routines:
  GetUnit

USE GridStatistics, ONLY : &
  !imported routines:
  GetMean

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

IMPLICIT NONE
! Global (i.e. public) Declarations:
! Global TYPE Definitions:

TYPE Aquifer
    TYPE(grid_real)   :: top
	TYPE(grid_real)   :: bottom
	TYPE(grid_real)   :: Ks !hydraulic conductivity (m/s)
	TYPE(grid_real)   :: sy !specific yield
	TYPE(grid_real)   :: KsAquitard !hydraulic conductivity of the aquitard
	TYPE(grid_integer):: domainBC !set domain and boundary conditions type 
	TYPE(grid_real)   :: valueBC ! value to be used for boundary condition  
    TYPE(grid_real)   :: head0 !head at time t-1
    TYPE(grid_real)   :: head1 !head at time t
END TYPE Aquifer

TYPE GroundwaterBasin
    INTEGER (KIND = short) :: nAquifers  !!number of aquifers in the basin
    TYPE (Aquifer), POINTER :: aquifer (:)
END TYPE GroundwaterBasin

!Global routines:
PUBLIC :: GroundwaterInit
PUBLIC :: GroundwaterUpdate
PUBLIC :: GroundwaterPointInit
PUBLIC :: GroundwaterOutputInit
PUBLIC :: GroundwaterOutput
PUBLIC :: GroundwaterRiverUpdate

!Global variables:
TYPE (GroundwaterBasin) :: basin
INTEGER (KIND = short)  :: dtGroundwater
TYPE (grid_real) :: vadoseDepth
TYPE (grid_real) :: groundwaterToRiver  !!discharge from groundwater to river (m3/s)
TYPE (grid_real) :: riverToGroundwater  !!discharge from river to groundwater (m3/s)
LOGICAL :: riverGroundwaterInteract

!Local (private) parameters:
INTEGER (KIND = short), PARAMETER, PRIVATE :: BC_DIRICHLET  = 2  ! head boundary condition
INTEGER (KIND = short), PARAMETER, PRIVATE :: BC_NEUMANN    = 3  ! flux boundary condition
INTEGER (KIND = short), PARAMETER, PRIVATE :: ACTIVE_CELL   = 1  ! active cell within the groundwater domain

!local variables
TYPE (grid_real), PRIVATE :: transmissivity  !!local transmissivity (m2/s)
TYPE (grid_real), PRIVATE :: fluxUpward  !!flux from lower aquitard (m/s)
TYPE (grid_real), PRIVATE :: fluxDownward  !!flux from upper aquitard (m/s)
TYPE (grid_real), PRIVATE :: riverDem !!reference digital elevation model to compute
                                        !!river water surface elevation (m asl)
TYPE (grid_real), PRIVATE :: streambedThickness !! streambed thickness (m)
TYPE (grid_real), PRIVATE :: streambedConductivity !! streambed conductivity (m/s)
TYPE (grid_integer), PRIVATE :: riverGroundwaterID !!id of river reaches that interact with groundwater
TYPE (Table), PRIVATE :: riverGroundwaterParameters

TYPE (ObservationalNetwork),PRIVATE :: sites
INTEGER (KIND = short), ALLOCATABLE, PRIVATE :: fileUnitPointGW (:)
INTEGER (KIND = short), PRIVATE :: fileUnitOutGW
TYPE (DateTime), PRIVATE :: timePointExport
REAL (KIND = float), PRIVATE :: volumeRecharge !! recharge volume (m3)
REAL (KIND = float), PRIVATE :: volumeResidual !! balance error (m3)
REAL (KIND = float), PRIVATE :: volumeNeumann !!volume from Neumann boundary condition (m3)
REAL (KIND = float), PRIVATE :: volumeDirichlet !!volume from Dirichlet boundary condition (m3)
REAL (KIND = float), PRIVATE :: volumeChange !! change of stored water volume in groundwater (m3)
REAL (KIND = float), PRIVATE :: volumeRiverToGroundwater !! volume from river to groundwater (m3) 
REAL (KIND = float), PRIVATE :: volumeGroundwaterToRiver !! volume from groundwater to river (m3)



!local (private) routines:
PRIVATE :: GroundwaterParameterUpdate
PRIVATE :: GroundwaterPointExport
PRIVATE :: GroundwaterRiverInit


!=======
CONTAINS
!=======
! Define procedures contained in this module. 
  
!==============================================================================
!| Description: 
!  Initialize groundwater
SUBROUTINE GroundwaterInit   & 
  !
  (inifile)       

IMPLICIT NONE


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

!local declarations:
TYPE (IniList)         :: iniDB
INTEGER (KIND = short) :: k , i, j
CHARACTER (LEN = 100)  :: string
REAL (KIND = float)    :: scalar

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

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

!read the number of aquifers
IF ( KeyIsPresent ('aquifers', iniDB ) ) THEN
    basin % naquifers = IniReadInt ( 'aquifers', iniDB  ) 
ELSE
      CALL Catch ('error', 'Groundwater',   &
			      'aquifers not found in configuration file' )
END IF

!allocate aquifers
ALLOCATE ( basin % aquifer ( basin % naquifers ) )

!load parameters
DO k = 1, basin % naquifers
    string = 'aquifer_' // ToString (k)
    IF (SectionIsPresent(TRIM(string), iniDB)) THEN	
        
        !load domain and boundary condition id
        IF (SubSectionisPresent (subsection = 'boundary-condition-id', &
                            section = TRIM(string), iniDB = iniDB) ) THEN
            
            CALL GridByIni ( ini = iniDB, &
                             grid = basin % aquifer (k) % domainBC, &
                             section = TRIM(string), &
                             subsection = 'boundary-condition-id')
            
            
        ELSE
            CALL Catch ('error', 'Groundwater',   &
                        'missing boundary-condition-id in configuration file: ', &
                         argument = TRIM(string) )
        END IF
        
        
        !load boundary condition value

        IF (SubSectionisPresent (subsection = 'boundary-condition-value', &
                            section = TRIM(string), iniDB = iniDB) ) THEN
            IF (KeyIsPresent ('scalar', iniDB, TRIM(string), &
                                'boundary-condition-value' ) ) THEN
                
                    scalar = IniReadReal ('scalar', iniDB, &
                                    TRIM(string), 'boundary-condition-value' )
                
                    CALL NewGrid (basin % aquifer (k) % valueBC, &
                                    basin % aquifer (k) % domainBC, scalar)
            ELSE
            
                    CALL GridByIni ( ini = iniDB, &
                                grid = basin % aquifer (k) % valueBC, &
                                section = TRIM(string), &
                                subsection = 'boundary-condition-value')
            END IF
        ELSE
            CALL Catch ('error', 'Groundwater',   &
                        'missing boundary-condition-value in configuration file: ', &
                            argument = TRIM(string) )
        END IF
        
        
        !load initial head
        IF (SubSectionisPresent (subsection = 'initial-head', &
                            section = TRIM(string), iniDB = iniDB) ) THEN
            
            IF (KeyIsPresent ('scalar', iniDB, TRIM(string), &
                                 'initial-head' ) ) THEN
                
                scalar = IniReadReal ('scalar', iniDB, &
                                      TRIM(string), 'initial-head' )
                
                CALL NewGrid (basin % aquifer (k) % head0, &
                              basin % aquifer (k) % domainBC, scalar)
            ELSE
                CALL GridByIni ( ini = iniDB, &
                                grid = basin % aquifer (k) % head0, &
                                section = TRIM(string), &
                                subsection = 'initial-head')
            END IF
            
            CALL NewGrid (basin % aquifer (k) % head1, &
                              basin % aquifer (k) % domainBC, 0.)
           
            basin % aquifer (k) % head1 = basin % aquifer (k) % head0
            
            !head boundary condition overlay
            DO i = 1, basin % aquifer (k) % domainBC % idim
              DO j = 1, basin % aquifer (k) % domainBC % jdim
                 IF ( basin % aquifer (k) % domainBC % mat (i,j) == &
                     BC_DIRICHLET ) THEN
            
                     basin % aquifer (k) % head0 % mat (i,j) = &
                     basin % aquifer (k) % valueBC % mat (i,j)
                     
                     basin % aquifer (k) % head1 % mat (i,j) = &
                     basin % aquifer (k) % valueBC % mat (i,j)
                     
                     
                END IF
              END DO
            END DO
            

        ELSE
            CALL Catch ('error', 'Groundwater',   &
                        'missing initial-condition-head in configuration file: ', &
                            argument = TRIM(string) )
        END IF
        
        
        !load top elevation
        IF (SubSectionisPresent (subsection = 'top-elevation', &
                            section = TRIM(string), iniDB = iniDB) ) THEN
            IF (KeyIsPresent ('scalar', iniDB, TRIM(string), &
                              'top-elevation' ) ) THEN
                
                scalar = IniReadReal ('scalar', iniDB, &
                                      TRIM(string), 'top-elevation' )
                
                CALL NewGrid (basin % aquifer (k) % top, &
                              basin % aquifer (k) % domainBC, scalar)
                
            ELSE
            
                CALL GridByIni ( ini = iniDB, &
                                grid = basin % aquifer (k) % top, &
                                section = TRIM(string), &
                                subsection = 'top-elevation')
           END IF

        ELSE
            CALL Catch ('error', 'Groundwater',   &
                        'missing top-elevation in configuration file: ', &
                            argument = TRIM(string) )
        END IF
        
        
        !load bottom elevation
        IF (SubSectionisPresent (subsection = 'bottom-elevation', &
                            section = TRIM(string), iniDB = iniDB) ) THEN
            IF (KeyIsPresent ('scalar', iniDB, TRIM(string), &
                              'bottom-elevation' ) ) THEN
                
                scalar = IniReadReal ('scalar', iniDB, &
                                      TRIM(string), 'bottom-elevation' )
                
                CALL NewGrid (basin % aquifer (k) % bottom, &
                              basin % aquifer (k) % domainBC, scalar)
                
            ELSE
            
                CALL GridByIni ( ini = iniDB, &
                                grid = basin % aquifer (k) % bottom, &
                                section = TRIM(string), &
                                subsection = 'bottom-elevation')
           END IF

        ELSE
            CALL Catch ('error', 'Groundwater',   &
                        'missing bottom-elevation in configuration file: ', &
                            argument = TRIM(string) )
        END IF
        
        
        
        !load hydraulic conductivity
        IF (SubSectionisPresent (subsection = 'hydraulic-conductivity', &
                            section = TRIM(string), iniDB = iniDB) ) THEN
            IF (KeyIsPresent ('scalar', iniDB, TRIM(string), &
                              'hydraulic-conductivity' ) ) THEN
                
                scalar = IniReadReal ('scalar', iniDB, &
                                      TRIM(string), 'hydraulic-conductivity' )
                
                CALL NewGrid (basin % aquifer (k) % Ks, &
                              basin % aquifer (k) % domainBC, scalar)
                
            ELSE
            
                CALL GridByIni ( ini = iniDB, &
                                grid = basin % aquifer (k) % Ks, &
                                section = TRIM(string), &
                                subsection = 'hydraulic-conductivity')
           END IF

        ELSE
            CALL Catch ('error', 'Groundwater',   &
                        'missing hydraulic-conductivity in configuration file: ', &
                            argument = TRIM(string) )
        END IF
        
        
        !load specific yield
        IF (SubSectionisPresent (subsection = 'specific-yield', &
                            section = TRIM(string), iniDB = iniDB) ) THEN
            IF (KeyIsPresent ('scalar', iniDB, TRIM(string), &
                              'specific-yield' ) ) THEN
                
                scalar = IniReadReal ('scalar', iniDB, &
                                      TRIM(string), 'specific-yield' )
                
                CALL NewGrid (basin % aquifer (k) % sy, &
                              basin % aquifer (k) % domainBC, scalar)
                
            ELSE
            
                CALL GridByIni ( ini = iniDB, &
                                grid = basin % aquifer (k) % sy, &
                                section = TRIM(string), &
                                subsection = 'specific-yield')
           END IF

        ELSE
            CALL Catch ('error', 'Groundwater',   &
                        'missing specific-yield in configuration file: ', &
                            argument = TRIM(string) )
        END IF
        
        
        ! aquitard hydraulic conductivity
        IF ( k <  basin % naquifers ) THEN 
            IF (SubSectionisPresent (subsection = 'aquitard-conductivity', &
                                section = TRIM(string), iniDB = iniDB) ) THEN
                IF (KeyIsPresent ('scalar', iniDB, TRIM(string), &
                                  'aquitard-conductivity' ) ) THEN
                
                    scalar = IniReadReal ('scalar', iniDB, &
                                          TRIM(string), 'aquitard-conductivity' )
                
                    CALL NewGrid (basin % aquifer (k) % KsAquitard, &
                                  basin % aquifer (k) % domainBC, scalar)
                
                ELSE
            
                    CALL GridByIni ( ini = iniDB, &
                                    grid = basin % aquifer (k) % KsAquitard, &
                                    section = TRIM(string), &
                                    subsection = 'aquitard-conductivity')
               END IF

            ELSE
                CALL Catch ('error', 'Groundwater',   &
                            'missing aquitard-conductivity in configuration file: ', &
                                argument = TRIM(string) )
            END IF
        END IF

        
        
    ELSE
        CALL Catch ('error', 'Groundwater',   &
			        'missing section in configuration file: ', &
                    argument = TRIM(string) )
    END IF 
    
END DO

!allocate variables
CALL NewGrid ( transmissivity, basin % aquifer (1) % domainBC, 0. )
CALL NewGrid ( fluxUpward, basin % aquifer (1) % domainBC, 0. )
CALL NewGrid ( fluxDownward, basin % aquifer (1) % domainBC, 0. )
CALL NewGrid ( vadoseDepth, basin % aquifer (1) % domainBC, 0. )

!check river-groundwater interaction
IF ( SectionIsPresent ( 'river-groundwater', iniDB) ) THEN
    riverGroundwaterInteract = .TRUE.
    CALL GroundwaterRiverInit ( inifile )
ELSE
    riverGroundwaterInteract = .FALSE.
END IF


!destroy iniDB
CALL IniClose (IniDB)


  END SUBROUTINE GroundwaterInit
  
  
  
  
!==============================================================================
!| Description: 
!  update groundwater head with Darcy law and mass conservation
!  in two dimensions, with a macrocopic cellular automata approach
!```
!            +-----+
!            |  N  |
!      +-----+-----+-----+
!      |  W  |  C  |  E  |
!      +-----+-----+-----+
!            |  S  |
!            +-----+
!```
! References:
!   Ravazzani, G., Rametta, D., & Mancini, M.. (2011). Macroscopic cellular 
!   automata for groundwater modelling: a first approach. Environmental 
!   modelling & software, 26(5), 634–643.
!
SUBROUTINE GroundwaterUpdate   & 
  !
  ( time, hillslopeFlux, percolation, capflux )       

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime),  INTENT(IN) :: time
TYPE (grid_real), INTENT(IN) :: hillslopeFlux !!flux from hillslope (m3/s)
TYPE (grid_real), INTENT(IN) :: percolation !! depp infiltration from soil balance (m/s)
TYPE (grid_real), INTENT(IN) :: capflux !! capillary flux from soil balance (m/s)


!local declarations: 
INTEGER (KIND = short) :: i, j, k
REAL (KIND = float) :: rechargeQ !!vertical recharge (m3/s)
REAL (KIND = float) :: transNC !!mean transmissivity North-Centre (m2/s)
REAL (KIND = float) :: transEC !!mean transmissivity East-Centre (m2/s)
REAL (KIND = float) :: transSC !!mean transmissivity South-Centre (m2/s)
REAL (KIND = float) :: transWC !!mean transmissivity West-Centre (m2/s)
REAL (KIND = float) :: fluxNC  !!flux North-Centre (m3/s)
REAL (KIND = float) :: fluxEC  !!flux East-Centre (m3/s)
REAL (KIND = float) :: fluxSC  !!flux South-Centre (m3/s)
REAL (KIND = float) :: fluxWC  !!flux West-Centre (m3/s)
REAL (KIND = float) :: fluxNESW !! total flux from the 4 directions (m3/s)
REAL (KIND = float) :: areaOfCell !! cell area (m2)
REAL (KIND = float) :: rToG !!river to groundwater discharge (m3/s)
REAL (KIND = float) :: gToR !!groundwater to river discharge (m3/s)

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

!reset counter variables
volumeNeumann = 0.
volumeDirichlet = 0.
volumeRecharge = 0.
volumeChange = 0.
volumeResidual = 0.

!update boundary condition
CALL GroundwaterParameterUpdate (time)


!save head1 into head0
DO k = 1, basin % nAquifers
    basin % aquifer (k) % head0 =  basin % aquifer (k) % head1
END DO

!update flux from hillslope of first (surface) aquifer
DO i = 1, basin % aquifer (1) % domainBC % idim
    DO j = 1, basin % aquifer (1) % domainBC % jdim
        IF ( basin % aquifer (1) % domainBC % mat (i,j) == &
             BC_NEUMANN ) THEN
            
            basin % aquifer (1) % valueBC % mat (i,j) = &
              hillslopeFlux % mat (i,j)
            
        END IF
    END DO
END DO

!update aquifers head
DO k = 1, basin % nAquifers
    
    !compute local transmissivity
    DO i = 1, basin % aquifer (k) % domainBC % idim
      DO j = 1, basin % aquifer (k) % domainBC % jdim
          IF ( basin % aquifer (k) % domainBC % mat(i,j) /=  &
			     basin % aquifer (k) % domainBC % nodata ) THEN
              
             !transmissivity: for confined aquifer = ks * aquifer depth
             transmissivity % mat (i,j) = MIN ( &
                 basin % aquifer (k) % Ks % mat (i,j) * &
                ( basin % aquifer (k) % head0 % mat (i,j) - &
                  basin % aquifer (k) % bottom % mat (i,j) ) , &

                  basin % aquifer (k) % Ks % mat (i,j) * &
                ( basin % aquifer (k) % top % mat (i,j) - &
                  basin % aquifer (k) % bottom % mat (i,j) )   &
								               )
                 
                 IF ( transmissivity % mat (i,j) <= 0. ) THEN
                     transmissivity % mat (i,j) = 0.01
                 END IF
          END IF
      END DO
    END DO
    
    
    !update head1
    DO i = 1, basin % aquifer (k) % domainBC % idim
      DO j = 1, basin % aquifer (k) % domainBC % jdim
          IF ( basin % aquifer (k) % domainBC % mat(i,j) == &
               ACTIVE_CELL ) THEN
              
             !harmonic mean transmissivity
             transNC = 2.0 * transmissivity % mat(i-1,j) * &
                             transmissivity % mat(i,j) / &
			               ( transmissivity  %mat(i-1,j) + &
                             transmissivity  %mat(i,j) )
             transSC = 2.0 * transmissivity % mat(i+1,j) * &
                             transmissivity % mat(i,j) / &
			               ( transmissivity % mat(i+1,j) + &
                             transmissivity % mat(i,j) )
             transEC = 2.0 * transmissivity % mat(i,j+1) * &
                             transmissivity % mat(i,j) / &
			               ( transmissivity % mat(i,j+1) + &
                             transmissivity % mat(i,j) )
             transWC = 2.0 * transmissivity % mat(i,j-1) * &
                             transmissivity % mat(i,j) / &
			               ( transmissivity % mat(i,j-1) + &
                             transmissivity % mat(i,j) )
             
             !flux from North
             IF ( basin % aquifer (k) % domainBC % mat(i-1,j) == &
                  BC_NEUMANN ) THEN
                 fluxNC = basin % aquifer (k) % valueBC % mat(i-1,j)
                 volumeNeumann = volumeNeumann + fluxNC * dtGroundwater
             ELSE
                 fluxNC = transNC * &
                     ( basin % aquifer (k) % head0 % mat (i-1,j) - &
                       basin % aquifer (k) % head0 % mat (i,j) )
                 
                 IF ( basin % aquifer (k) % domainBC % mat(i-1,j) == &
                      BC_DIRICHLET ) THEN
                     volumeDirichlet = volumeDirichlet + &
                                       fluxNC * dtGroundwater
                 END IF
             END IF
             
             !flux from East
             IF ( basin % aquifer (k) % domainBC % mat(i,j+1) == &
                  BC_NEUMANN ) THEN
                 fluxEC = basin % aquifer (k) % valueBC % mat(i,j+1)
                 volumeNeumann = volumeNeumann + fluxEC * dtGroundwater
             ELSE
                 fluxEC = transEC * &
                     ( basin % aquifer (k) % head0 % mat (i,j+1) - &
                       basin % aquifer (k) % head0 % mat (i,j) )
                 
                 IF ( basin % aquifer (k) % domainBC % mat(i,j+1) == &
                      BC_DIRICHLET ) THEN
                     volumeDirichlet = volumeDirichlet + &
                                       fluxEC * dtGroundwater
                 END IF
             END IF
             
             !flux from South
             IF ( basin % aquifer (k) % domainBC % mat(i+1,j) == &
                  BC_NEUMANN ) THEN
                 fluxSC = basin % aquifer (k) % valueBC % mat(i+1,j)
                 volumeNeumann = volumeNeumann + fluxSC * dtGroundwater
             ELSE
                 fluxSC = transSC * &
                     ( basin % aquifer (k) % head0 % mat (i+1,j) - &
                       basin % aquifer (k) % head0 % mat (i,j) )
                 
                 IF ( basin % aquifer (k) % domainBC % mat(i+1,j) == &
                      BC_DIRICHLET ) THEN
                     volumeDirichlet = volumeDirichlet + &
                                       fluxSC * dtGroundwater
                 END IF
             END IF
             
             !flux from West
             IF ( basin % aquifer (k) % domainBC % mat(i,j-1) == &
                  BC_NEUMANN ) THEN
                 fluxWC = basin % aquifer (k) % valueBC % mat(i,j-1)
                 volumeNeumann = volumeNeumann + fluxWC * dtGroundwater
             ELSE
                 fluxWC = transWC * &
                     ( basin % aquifer (k) % head0 % mat (i,j-1) - &
                       basin % aquifer (k) % head0 % mat (i,j) )
                 
                 IF ( basin % aquifer (k) % domainBC % mat(i,j-1) == &
                      BC_DIRICHLET ) THEN
                     volumeDirichlet = volumeDirichlet + &
                                       fluxWC * dtGroundwater
                 END IF
             END IF
             
             !total flux
             fluxNESW = fluxNC + fluxEC + fluxSC + fluxWC
             
             !area of the cell
             areaOfCell = CellArea (transmissivity, i, j) 
             
             !vertical recharge
             IF ( ALLOCATED (percolation % mat ) ) THEN
                 IF ( k == 1 .AND. percolation % mat (i,j) /= &
                                   percolation % nodata ) THEN
                     rechargeQ = ( percolation % mat (i,j) - &
                                   capflux % mat (i,j) ) * areaOfCell
                     volumeRecharge = volumeRecharge + &
                                      rechargeQ * dtGroundwater
                 ELSE
                     rechargeQ = 0.
                 END IF
             ELSE
                 rechargeQ = 0.
             END IF
             
             !upward flux from lower aquifer
             IF (k /= basin % nAquifers ) THEN
			    fluxUpward % mat(i,j) = &
                    basin % aquifer (k) % ksAquitard % mat(i,j) * &
                ( ( basin % aquifer (k+1) % head0 % mat(i,j) -  &
                    basin % aquifer (k) % head0 % mat(i,j) ) / &
                  ( basin % aquifer (k) % bottom % mat(i,j) - &
                    basin % aquifer (k+1) % top % mat(i,j)  ) )
             ELSE
			    fluxUpward % mat (i,j) = 0.0
             END IF
             
             !update head1  
             IF ( k == 1) THEN
                 fluxDownward % mat (i,j) = 0.
                 IF ( riverGroundwaterInteract ) THEN
                     rToG = riverToGroundwater % mat (i,j)
                     gToR = - groundwaterToRiver % mat (i,j)
                 ELSE
                     rToG = 0.
                     gToR = 0.
                 END IF
             ELSE
                 rToG = 0.
                 gToR = 0.
             END IF
             
             basin % aquifer (k) % head1 % mat(i,j) = &
                 basin % aquifer (k) % head0 % mat(i,j) + &
                 ( fluxNESW + rechargeQ + rToG + gToR ) * &
                 dtGroundwater / areaOfCell / &
                 basin % aquifer (k) % sy % mat(i,j) + &
                 ( fluxUpward % mat(i,j) + fluxDownward % mat(i,j) ) * &
                 dtGroundwater / basin % aquifer (k) % sy % mat(i,j)
             
             IF ( ISNAN ( basin % aquifer (k) % head1 % mat(i,j) ) )  THEN
               basin % aquifer (k) % head1 % mat(i,j) = &
                   basin % aquifer (k) % head0 % mat(i,j)
                volumeChange = 0.
             ELSE
                !volume change
                 volumeChange = volumeChange + &
                 ( basin % aquifer (k) % head1 % mat(i,j) - &
                   basin % aquifer (k) % head0 % mat(i,j) ) * &
                   basin % aquifer (k) % sy % mat(i,j) * areaOfCell 
             END IF
             
             !swap fluxUpward and fluxDownward
             fluxDownward % mat (i,j) = - fluxUpward % mat (i,j) 
             
             !update vadose zone depth
             IF ( k == 1) THEN
                 vadoseDepth % mat (i,j) = & 
                     basin % aquifer (k) % top % mat(i,j) - &
                     basin % aquifer (k) % head1 % mat(i,j)
             END IF
             
          END IF
      END DO
    END DO
    
END DO

!compute balance error
volumeResidual = volumeRecharge + volumeDirichlet + &
                 volumeNeumann  + volumeRiverToGroundwater - &
                 volumeGroundwaterToRiver - volumeChange

!write point data
IF ( time == timePointExport ) THEN
   CALL GroundwaterPointExport ( time )
END IF


RETURN
END SUBROUTINE GroundwaterUpdate


    
    
  
!==============================================================================
!| Description: 
!  update boundary condition map that change in time
SUBROUTINE GroundwaterParameterUpdate   & 
  !
  (time)       

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT (IN) :: time

!local declarations:
CHARACTER (LEN = 300) :: filename
CHARACTER (LEN = 300) :: varname
INTEGER (KIND = short) :: k, i, j

!------------------------------end of declarations-----------------------------
  
  
  
!boundary condition
DO k = 1, basin % naquifers
    IF (  time == basin % aquifer (k) % valueBC % next_time ) THEN
       !destroy current grid
       filename = basin % aquifer (k) % valueBC % file_name
       varname = basin % aquifer (k) % valueBC % var_name
       CALL GridDestroy (basin % aquifer (k) % valueBC )
       !read new grid in netcdf file
       CALL NewGrid (basin % aquifer (k) % valueBC, TRIM(filename), NET_CDF, &
                      variable = TRIM(varname), time = time)
       
       !head boundary condition overlay
        DO i = 1, basin % aquifer (k) % domainBC % idim
            DO j = 1, basin % aquifer (k) % domainBC % jdim
                IF ( basin % aquifer (k) % domainBC % mat (i,j) == &
                    BC_DIRICHLET ) THEN
            
                    basin % aquifer (k) % head0 % mat (i,j) = &
                    basin % aquifer (k) % valueBC % mat (i,j)
                     
                    basin % aquifer (k) % head1 % mat (i,j) = &
                    basin % aquifer (k) % valueBC % mat (i,j)
               END IF
            END DO
        END DO
    END IF
END DO
  
RETURN
END SUBROUTINE GroundwaterParameterUpdate  
    
    
!==============================================================================
!| Description:
!   Initialize export of point site data
SUBROUTINE GroundwaterPointInit &
!
( pointfile, path_out, time )

IMPLICIT NONE

!Arguments with intent (in):
CHARACTER (LEN = *), INTENT(IN) :: pointfile  !!file containing coordinate of points
CHARACTER (LEN = *), INTENT(IN) :: path_out  !!path of output folder
TYPE (DateTime),     INTENT(IN) :: time  !!start time

!local declarations
INTEGER (KIND = short) :: err_io
INTEGER (KIND = short) :: fileunit
INTEGER (KIND = short) :: k

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

timePointExport = time

fileunit = GetUnit ()
OPEN ( unit = fileunit, file = pointfile(1:LEN_TRIM(pointfile)), &
      status='OLD', iostat = err_io )

IF ( err_io /= 0 ) THEN
	!file does not exist
    CALL Catch ('error', 'Groundwater', 'out point file does not exist')
END IF 

!Read metadata
CALL ReadMetadata (fileunit, sites)

!check dt
IF (.NOT. ( MOD ( sites % timeIncrement, dtGroundwater ) == 0 ) ) THEN
  CALL Catch ('error', 'Groundwater', 'dt out sites must be multiple of dtGroundwater ')
END IF

CLOSE ( fileunit )

sites % description = 'groundwater head data exported from FEST'

sites % unit = 'm'

sites % offsetZ = 0.

!allocate file unit
ALLOCATE ( fileUnitPointGW ( basin % nAquifers ) )

!open and initialize files
DO k = 1, basin % nAquifers
    fileUnitPointGW (k) = GetUnit ()
    OPEN ( unit = fileUnitPointGW (k), &
    file = TRIM(path_out) // 'point_aquifer_' //  TRIM(ToString (k)) // '.fts' )
    
    CALL WriteMetadata ( network = sites, fileunit = fileUnitPointGW (k) )

    CALL WriteData (sites, fileUnitPointGW (k), .TRUE.)  
    
END DO 

RETURN
END SUBROUTINE GroundwaterPointInit    


!==============================================================================
!| Description:
!   Export of point site data
SUBROUTINE GroundwaterPointExport &
!
( time )

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT (IN) :: time

!local declarations:
INTEGER (KIND = short) :: k

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

!set current time
sites % time = time 

DO k = 1, basin % nAquifers
    
    !populate data
    CALL AssignDataFromGrid ( basin % aquifer (k) % head1, sites )
    
    !write data
    CALL WriteData ( sites, fileUnitPointGW (k) )
    
END DO
    
timePointExport = timePointExport + sites % timeIncrement


RETURN
END SUBROUTINE GroundwaterPointExport



!==============================================================================
!| Description:
!   Initialize files for exporting output results
SUBROUTINE GroundwaterOutputInit &
!
( path_out )

IMPLICIT NONE

!Arguments with intent (in):
CHARACTER (LEN = *), INTENT(IN) :: path_out  !!path of output folder

!local declarations
REAL (KIND = float) :: area !! (km2)
INTEGER (KIND = short) :: i, j
INTEGER (KIND = short) :: nvars



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

!compute groundwater extent area
area = 0.
DO i = 1, basin % aquifer (1) % domainBC % idim
    DO j = 1, basin % aquifer (1) % domainBC % jdim
        IF (basin % aquifer (1) % domainBC % mat (i,j) == ACTIVE_CELL ) THEN
            area =  area + CellArea ( basin % aquifer (1) % domainBC, i, j )
        END IF
    END DO
END DO
area = area / 1000000. !conversion m2 -> km2

!number of variables to write in output file
nvars = basin % nAquifers + 5

IF ( riverGroundwaterInteract ) THEN
    nvars = nvars + 2
END IF


!open and initialize file
fileUnitOutGW = GetUnit ()
OPEN ( unit = fileUnitOutGW, file = TRIM(path_out) // 'aquifer.out' )
WRITE (fileUnitOutGW,'(a)') 'Output variables of groundwater simulation'
WRITE (fileUnitOutGW,fmt = '(a, F12.2)') 'extent area (km2): ', area
WRITE (fileUnitOutGW,fmt = '(a, I3)') 'number of variables: ', nvars
WRITE (fileUnitOutGW,*)
WRITE (fileUnitOutGW,*)
WRITE (fileUnitOutGW,'(a)') 'data'
WRITE (fileUnitOutGW,'(a)', advance='no') 'DateTime     '

DO i = 1, basin % nAquifers
    WRITE (fileUnitOutGW,'(a)', advance='no') 'head_aquifer_' // &
                                          TRIM(ToString(i)) // '_m'
END DO

IF ( riverGroundwaterInteract ) THEN

    WRITE (fileUnitOutGW,'(a)') '   recharge_m3   volume_neumann_m3   &
           volume_dirichlet_m3    volume_change_m3    residual_m3  & 
           river_to_groundwater_m3    groundwater_to_river_m3 '

ELSE
    WRITE (fileUnitOutGW,'(a)') '   recharge_m3   volume_neumann_m3   &
           volume_dirichlet_m3    volume_change_m3    residual_m3 '
END IF
RETURN
END SUBROUTINE GroundwaterOutputInit



!==============================================================================
!! Description:
!!   Write results on output file
SUBROUTINE GroundwaterOutput &
!
( time )

IMPLICIT NONE

!Arguments with intent (in):
TYPE (DateTime), INTENT(IN) :: time  !! current time

!local variables:
INTEGER (KIND = short) :: k
REAL (KIND = float) :: meanHead

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

timeString = time
WRITE (fileUnitOutGW,'(a)',advance='no') timeString

DO k = 1, basin % nAquifers
    meanHead = GetMean ( basin % aquifer (k) % head1,  &
                         maskInteger = basin % aquifer (k) % domainBC ) 
    WRITE (fileUnitOutGW,'(6x,E12.5)',advance='no') meanHead
END DO

IF ( riverGroundwaterInteract ) THEN
    WRITE (fileUnitOutGW,'(7(5x,E12.5))') volumeRecharge, volumeNeumann, &
                            volumeDirichlet, volumeChange, volumeResidual, &
                            volumeRiverToGroundwater, volumeGroundwaterToRiver

ELSE

    WRITE (fileUnitOutGW,'(5(5x,E12.5))') volumeRecharge, volumeNeumann, &
                            volumeDirichlet, volumeChange, volumeResidual 
END IF

RETURN
END SUBROUTINE GroundwaterOutput

!==============================================================================
!| Description:
!   Configure river-groundwater interaction
SUBROUTINE GroundwaterRiverInit &
!
( inifile )

IMPLICIT NONE

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

!local variables:
TYPE (IniList)         :: iniDB
INTEGER (KIND = short) :: i,j


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

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

!load river id
IF (SubSectionisPresent (subsection = 'river-id', &
                            section = 'river-groundwater', iniDB = iniDB) ) THEN
     CALL GridByIni ( ini = iniDB, &
                             grid = riverGroundwaterID, &
                             section = 'river-groundwater', &
                             subsection = 'river-id')
ELSE
     CALL Catch ('error', 'Groundwater',   &
                'missing river-id in configuration file' )
END IF

!load river dem
IF (SubSectionisPresent (subsection = 'river-dem', &
                            section = 'river-groundwater', iniDB = iniDB) ) THEN
     CALL GridByIni ( ini = iniDB, &
                             grid = riverDem, &
                             section = 'river-groundwater', &
                             subsection = 'river-dem')
ELSE
     CALL Catch ('error', 'Groundwater',   &
                'missing river-dem in configuration file' )
END IF

!exchange parameters
CALL TableNew ( file = inifile, tab = riverGroundwaterParameters, &
                id = 'river-groundwater')

!allocate exchange flux maps
CALL NewGrid ( riverToGroundwater, riverGroundwaterID, 0. )
CALL NewGrid ( groundwaterToRiver, riverGroundwaterID, 0. )

!allocate streambed parameter maps
CALL NewGrid ( streambedThickness, riverGroundwaterID, 0. )
CALL NewGrid ( streambedConductivity, riverGroundwaterID, 0. )

!populate streambed parameter maps
DO i = 1, riverGroundwaterID % idim
    DO j = 1, riverGroundwaterID % jdim
        IF ( riverGroundwaterID % mat (i,j) /= &
             riverGroundwaterID % nodata ) THEN
            
            CALL TableGetValue ( &
                valueIn =  REAL(riverGroundwaterID % mat (i,j)),&
                tab = riverGroundwaterParameters, &
                keyIn = 'id', keyOut ='streambed-conductivity', &
                match = 'exact', valueOut = streambedConductivity % mat (i,j) )
            
            CALL TableGetValue ( &
                valueIn =  REAL(riverGroundwaterID % mat (i,j)),&
                tab = riverGroundwaterParameters, &
                keyIn = 'id', keyOut ='streambed-thickness', &
                match = 'exact', valueOut = streambedThickness % mat (i,j) )
            
        END IF
    END DO
END DO


!free db
CALL IniClose ( iniDB )


RETURN
END SUBROUTINE GroundwaterRiverInit



!==============================================================================
!| Description:
!   Update river-groundwater exchange fluxes
SUBROUTINE GroundwaterRiverUpdate &
!
( waterDepth, topWidth )

IMPLICIT NONE

!Arguments with intent (in):
TYPE (grid_real), INTENT(IN) :: waterDepth !!river water depth (m)
TYPE (grid_real), INTENT(IN) :: topWidth !!river top width (m)


!local variables:
INTEGER (KIND = short) :: i, j
REAL (KIND = float) :: riverWSE !river water surface elevation (m asl)
!TYPE (grid_real) :: waterTable  !!aquifer water table (m)
REAL (KIND = float) :: waterTable
REAL (KIND = float) :: area

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

volumeRiverToGroundwater = 0.
volumeGroundwaterToRiver = 0.

!waterTable = basin % aquifer (1) % head1

DO i = 1, riverGroundwaterID % idim
    DO j = 1, riverGroundwaterID % jdim
        IF ( riverGroundwaterID % mat (i,j) /= &
             riverGroundwaterID % nodata ) THEN 
            waterTable = basin % aquifer (1) % head1 % mat (i,j)
            riverWSE = waterDepth % mat (i,j) + riverDem % mat (i,j)
            area = CellArea (riverGroundwaterID, i, j)
            
            IF ( waterTable > riverWSE ) THEN
                
                groundwaterToRiver % mat (i,j) = &
                   ( waterTable - riverWSE )              * & !head difference
                     streambedConductivity % mat (i,j)    / & !conductivity
                     streambedThickness % mat (i,j)       * & !thickness
                     topWidth % mat (i,j)                 * & !river width
                     area ** 0.5                             !cell size
                
                riverToGroundwater % mat (i,j) = 0.
                
            ELSE  IF ( waterTable  <  riverWSE .AND. &
	                   waterTable  >  riverDem % mat (i,j) ) THEN
                
               ! when waterdepth > 10 cm, compute river discharge toward groundwater  
	           IF ( waterDepth % mat (i,j) > 0.10 ) THEN   
                   riverToGroundwater % mat(i,j) =             &
                       ( riverWSE - waterTable )             * & !head difference
                         streambedConductivity % mat (i,j)   / & !conductivity
                         streambedThickness % mat (i,j)      * & !thickness
                         topWidth % mat (i,j)                * & !river width
                         area ** 0.5                             !cell size
                   
                   groundwaterToRiver % mat (i,j) = 0.
			   ELSE
			       riverToGroundwater % mat(i,j) = 0.
               END IF
               
          ELSE  IF ( waterTable < riverDem % mat (i,j) ) THEN 
	            ! when waterdepth > 10 cm, compute river discharge toward groundwater   
	           IF (  waterDepth % mat (i,j) > 0.10 ) THEN  
	                riverToGroundwater % mat(i,j) =            &      
                         waterDepth % mat (i,j)              * & !head
                         streambedConductivity % mat (i,j)   / & !conductivity
                         streambedThickness % mat (i,j)      * & !thickness
                         topWidth % mat (i,j)                * & !river width
                         area ** 0.5      !cell size
                    
                    groundwaterToRiver % mat (i,j) = 0.
        
			   ELSE
			       riverToGroundwater % mat(i,j) = 0.
               END IF
                
            
          END IF
            
          volumeRiverToGroundwater = volumeRiverToGroundwater + &
                                     riverToGroundwater % mat (i,j) * &
                                     dtGroundwater
          volumeGroundwaterToRiver = volumeGroundwaterToRiver + &
                                     groundwaterToRiver % mat (i,j) * &
                                     dtGroundwater
            
            
        END IF
    END DO
END DO


RETURN
END SUBROUTINE GroundwaterRiverUpdate
    
END MODULE Groundwater