!! Solve soil water and energy balance 
!|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.6 - 20th July 2023  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 28/Jun/2010 | Initial code. chiara corbari and davide rabuffetti |
! | 1.1      | 03/May/2011 | Correct energy fluxes when snow is present |
! | 1.2      | 10/Nov/2016 | Module rewritten to add new infiltration models and different soil properties management |
! | 1.3      | 25/Jul/2022 | subsurface flow changed from muskingum to darcy. Runon option removed |
! | 1.4      | 27/Jan/2023 | bug correction: percolationcell set to zero on hillslope and rivers in SetSoilDepth |
! | 1.5      | 10/Feb/2023 | modified to not use river network |
! | 1.6      | 20/Jul/2023 | SetSoilDepth renamed to PercolationAndCaprise and modified |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
!   Solve soil water and energy balance to update soil moisture and
!   compute actual evapotranspiration, runoff, and groundwater 
!   recharge.
!   Soil column is assumed to be divided in two layers:
!
! * the surface layer where roots develop (root zone)
!
! * the lower unsaturated layer that transfers
!      water to recharge groundwater in landplain (transmission zone)
!
MODULE SoilBalance        

! Modules used: 

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

USE DomainProperties, ONLY : &
    !imported variables:
    mask

USE StringManipulation, ONLY: &
  !Imported routines:
  ToString

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

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

USE Chronos, ONLY: &
  !Imported definitions:
  DateTime, &
  !Imported operators and assignment:
  OPERATOR (>=), OPERATOR (<=), &
  ASSIGNMENT( = )

USE LogLib, ONLY: &
  !Imported routines:
  Catch

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

USE SoilProperties, ONLY : &
  !Imported definitions:
  !SoilType, &
  !imported routines:
  UnsHydCond, Psi

USE Infiltration, ONLY : &
    !imported routines:
    InfiltrationInit, &
    Sat2scn, SCScurveNumber, &
    Philip, GreenAmpt, &
    !imported variables:
    infiltrationModel, &
    soilDivisions, &
    curveNumber, &
    storativity, &
    sEff, ism, cuminf, &
    raincum, wiltingPoint, &
    fieldCapacity, thetar, &
    thetas, ksat, psdi, psic, &
    abstractionRatio, phy, &
    !imported parameters:
    SCS_CN, PHILIPEQ, &
    GREEN_AMPT, ROSS_BC, ROSS_VG 

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

USE Evapotranspiration, ONLY : &
    !imported routines
    ETinit, PETupdate, &
    !imported variables:
    pet, pe, pt

USE Morphology, ONLY: &
  !Imported routines:
  DownstreamCell, CheckOutlet

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

USE StringManipulation, ONLY: &
  !Imported routines:
  StringToShort

USE MorphologicalProperties, ONLY: &
  !imported variables:
  dem

USE Snow, ONLY : &
  !imported variables:
  swe, waterInSnow, dtSnow

USE Glacier, ONLY : &
  !imported variables:
  icewe, waterInIce, dtIce

!Global declarations:

!public routines
PUBLIC :: InitSoilBalance
PUBLIC :: SolveSoilBalance

!public variables:

INTEGER (KIND = short) :: dtSoilBalance

!TYPE(NetworkSoilBalance):: soilNetwork !!define spatial scheme to solve balance

!parameters
TYPE (grid_real)     :: ksat_sub ![m/s]
TYPE (grid_real)     :: soilDepth !!soil depth [m]
TYPE (grid_real)     :: soilDepthRZ !! root zone depth [m]
TYPE (grid_real)     :: soilDepthTZ !!transmission zone depth [m]

!common state variables
TYPE (grid_real)    :: soilSat !!mean soil relative saturation  [0-1]
TYPE (grid_real)    :: soilSatRZ !!root zone soil relative saturation  [0-1]
TYPE (grid_real)    :: soilSatTZ !!transmission zone soil relative saturation  [0-1]
TYPE (grid_real)    :: soilMoisture !!mean volumetric water content  [m3/m3]
TYPE (grid_real)    :: soilMoistureRZ !!root zone volumetric water content  [m3/m3]
TYPE (grid_real)    :: soilMoistureTZ !!transmission zone volumetric water content  [m3/m3]
TYPE (grid_integer) :: rainFlag !!set if it is raining in a cell 
TYPE (grid_integer) :: interstormDuration !!interstorm duration [s]
TYPE (grid_real)    :: QinSoilSub !!input discharge in soil subsurface [m3/s]
TYPE (grid_real)    :: QoutSoilSub !!output discharge from soil subsurface [m3/s]
TYPE (grid_real)    :: balanceError !!balance error [mm]
TYPE (grid_real)    :: deltaSoilMoisture !! time step soil mositure variation (m)


!vertical fluxes
TYPE (grid_real)    :: rainBalance  !!actual rainfall rate as input to soil balance [m/s]
TYPE (grid_real)    :: infilt !!infiltration rate [m/s]
TYPE (grid_real)    :: runoff !!runoff rate [m/s]
TYPE (grid_real)    :: et !! actual evapotranspiration rate [m/s]
TYPE (grid_real)    :: percolation !! deep percolation through soil [m/s]
TYPE (grid_real)    :: percolationFactor !! deep percolation factor [-]
TYPE (grid_real)    :: capRise !! capillary rise [m/s]

!energy balance fluxes
TYPE (grid_real)    :: Ts 
TYPE (grid_real)    :: Rnetta 
TYPE (grid_real)    :: Xle 
TYPE (grid_real)    :: Hse
TYPE (grid_real)    :: Ge 
TYPE (grid_real)    :: Ta_prec
TYPE (grid_real)    :: Ts_prec


!Local routines
PRIVATE :: WetlandIsWet
PRIVATE :: SetInitialCondition
PRIVATE :: SetWetland
PRIVATE :: PercolationAndCaprise
PRIVATE :: UpdateSoilMoisture

!local declarations:

!parameters
INTEGER, PRIVATE, PARAMETER  :: HILLSLOPE       = 0
INTEGER, PRIVATE, PARAMETER  :: CHANNEL         = 1
INTEGER, PRIVATE, PARAMETER  :: LAKE            = 2
INTEGER, PRIVATE, PARAMETER  :: LANDPLAIN       = 3


!wetland
TYPE (grid_integer), PRIVATE :: wetland
INTEGER, ALLOCATABLE, PRIVATE :: wetCode (:)
TYPE (DateTime), ALLOCATABLE, PRIVATE :: wetBegin (:)
TYPE (DateTime), ALLOCATABLE, PRIVATE :: wetEnd (:)

!balance
REAL (KIND = float), PRIVATE :: thresholdStartEvent !threshold to consider storm initiated [m/s]
REAL (KIND = float), PRIVATE :: interstorm !duration of interstorm period to terminate an event [s]
REAL (KIND = float), PRIVATE :: isd !!initial saturation degree, used for cold start
TYPE (grid_integer), PRIVATE :: balanceId  !!id code for solving water balance
LOGICAL            , PRIVATE :: saturatedByGroundwater !!groundwater table intercepts root zone
                
!=======
CONTAINS
!=======
! Define procedures contained in this module. 
  
!==============================================================================
!| Description: 
!  Initialize soil water balance 
SUBROUTINE InitSoilBalance   & 
  !
  (inifile, flowDirection, time)       

IMPLICIT NONE

!Arguments with intent(in):
CHARACTER (LEN = *), INTENT(IN) :: inifile !!stores configuration information
TYPE (grid_integer), INTENT(IN) :: flowDirection !!flow direction map 
TYPE (DateTime),     INTENT(IN) :: time !!start time

!local declarations:
TYPE (IniList) :: iniDB
CHARACTER (LEN = 1000) :: filename
INTEGER (KIND = short) :: k, i, j, iin, jin, is, js
REAL (KIND = float) :: scalar

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

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

!read balance id
IF (SectionIsPresent('balance-id', iniDB)) THEN
  CALL GridByIni (iniDB, balanceId, section = 'balance-id')
ELSE !grid is mandatory: stop the program if not present
   CALL Catch ('error', 'SoilBalance',   &
			   'error in loading balance-id: ' ,  &
			    argument = 'section not defined in ini file' )
END IF

!read subsurface soil hydraulic conductivity
IF (SectionIsPresent('ksat-subsurface', iniDB)) THEN
    IF (KeyIsPresent ('scalar', iniDB, 'ksat-subsurface') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'ksat-subsurface')
        CALL NewGrid (ksat_sub, mask, scalar)
    ELSE
        CALL GridByIni (iniDB, ksat_sub, section = 'ksat-subsurface')
    END IF
ELSE !grid is mandatory: stop the program if not present
   CALL Catch ('error', 'SoilBalance',   &
			   'error in loading subsurface conductivity: ' ,  &
			    argument = 'section not defined in ini file' )
END IF

!read soil depth (m)
IF (SectionIsPresent('soil-depth', iniDB)) THEN
    IF (KeyIsPresent ('scalar', iniDB, 'soil-depth') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'soil-depth')
        CALL NewGrid (soilDepth, mask, scalar)
    ELSE
        CALL GridByIni (iniDB, soilDepth, section = 'soil-depth')
    END IF
ELSE !grid is mandatory: stop the program if not present
   CALL Catch ('error', 'SoilBalance',   &
			   'error in loading soil depth: ' ,  &
			    argument = 'section not defined in ini file' )
END IF



!read deep percolation factor (-)
IF (SectionIsPresent('deep-percolation-factor', iniDB)) THEN
    IF (KeyIsPresent ('scalar', iniDB, 'deep-percolation-factor') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'deep-percolation-factor')
        CALL NewGrid (percolationFactor, mask, scalar)
    ELSE
        CALL GridByIni (iniDB, percolationFactor, section = 'deep-percolation-factor')
    END IF
ELSE !grid is optional: default = 1.
   CALL NewGrid (percolationFactor, mask, 1.0)
END IF



!read root zone depth (m)
IF (SectionIsPresent('root-zone-depth', iniDB)) THEN
    IF (KeyIsPresent ('scalar', iniDB, 'root-zone-depth') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'root-zone-depth')
        CALL NewGrid (soilDepthRZ, mask, scalar)
    ELSE
        CALL GridByIni (iniDB, soilDepthRZ, section = 'root-zone-depth')
    END IF
ELSE !grid is mandatory: stop the program if not present
   CALL Catch ('error', 'SoilBalance',   &
			   'error in loading root zone depth: ' ,  &
			    argument = 'section not defined in ini file' )
END IF


!compute transmission zone depth
CALL NewGrid ( soilDepthTZ, mask, 0.)

DO i = 1, mask % idim
    DO j = 1, mask % jdim
        IF ( mask % mat (i,j) /= mask % nodata ) THEN
            soilDepthTZ % mat(i,j) = soilDepth % mat(i,j) -  &
                                     soilDepthRZ % mat(i,j)
            IF ( soilDepthTZ % mat(i,j) <= 0. ) THEN
                soilDepthTZ % mat(i,j) = 0.1
            END IF
        END IF
    END DO
END DO

    
!configure evapotranspiration
IF (KeyIsPresent('evapotranspiration', iniDB)) THEN	
   filename = IniReadString ('evapotranspiration', iniDB)
   CALL ETinit (filename, time)
ELSE
    CALL Catch ('error', 'SoilBalance',   &
			      'evapotranspiration not found in configuration file' )
END IF


!threshold to initiate storm period, read and convert from mm/h to m/s
thresholdStartEvent = IniReadReal ('threshold-storm-start', iniDB)
thresholdStartEvent = thresholdStartEvent / 1000. / 3600.

!interstorm duration to terminate an event, read and convert from hours to seconds
interstorm = IniReadReal ('interstorm', iniDB)
interstorm = interstorm * 3600.
    

!configure infiltration
IF (KeyIsPresent('infiltration', iniDB)) THEN	
   filename = IniReadString ('infiltration', iniDB)
   CALL InfiltrationInit (filename, soilMoisture, soilDepth )
ELSE
    CALL Catch ('error', 'SoilBalance',   &
			      'infiltration not found in configuration file' )
END IF

  
!set initial condition
CALL SetInitialCondition (iniDB)

  
!set wetland
CALL SetWetland (iniDB, inifile)
  
!allocate coomon variables used to solve soil water and energy balance
!vertical fluxes
CALL NewGrid (rainBalance, mask, 0.)
CALL NewGrid (infilt, mask, 0.)
CALL NewGrid (runoff, mask, 0.)
CALL NewGrid (et, mask, 0.)
CALL NewGrid (percolation, mask, 0.)
CALL NewGrid (capRise, mask, 0.)
          
!energy balance
CALL NewGrid (Ts, mask, 0.)
CALL NewGrid (Rnetta, mask, 0.)
CALL NewGrid (Xle, mask, 0.)
CALL NewGrid (Hse, mask, 0.)
CALL NewGrid (Ge, mask, 0.) 
CALL NewGrid (Ta_prec, mask, 0.) 
CALL NewGrid (Ts_prec, mask, 0.) 

!lateral fluxes at time t
CALL NewGrid (QinSoilSub, mask, 0.)
CALL NewGrid (QoutSoilSub, mask, 0.)

!balance error
CALL NewGrid (balanceError, mask, 0.)

!soil mositure variation
CALL NewGrid (deltaSoilMoisture, mask, 0.)


!  Configuration terminated. Deallocate ini database
CALL IniClose (iniDB)  
  
RETURN
END SUBROUTINE InitSoilBalance

  
!==============================================================================
!| Description: 
!  Solve soil water balance 
SUBROUTINE SolveSoilBalance   & 
  !
  (time, rain, irrigation, flowdir, vf, vadose)       

IMPLICIT NONE

!Arguments with intent(in):
!INTEGER (KIND = short), INTENT(IN) :: dt !!computation time step [s]
TYPE (DateTime), INTENT(IN) :: time !current date and time
TYPE (grid_real), INTENT(IN) :: rain !!rainfall rate [m/s]
TYPE (grid_real), INTENT(IN) :: irrigation !!irrigation rate [m/s]
TYPE (grid_integer), INTENT(IN) :: flowdir
TYPE (grid_real), INTENT(IN) :: vf !!vegetation fraction [0-1]
TYPE (grid_real), INTENT(IN) :: vadose !vadose zone depth [m]


!local declarations:
INTEGER (KIND = short) :: k, iin, jin, is, js, id, i, j
REAL (KIND = double) :: runoffcell !runoff of single cell [m/s]
REAL (KIND = float) :: percolationcellRZ !root zone percolation of single cell [m/s]
REAL (KIND = float) :: percolationcellTZ !transmission zone percolation of single cell [m/s]
REAL (KIND = float) :: Tdeep !deep surface temperature in lake
REAL (KIND = float) :: alpha, beta !!used to compute actual evapotranspiration [0-1]
REAL (KIND = float) :: sm, wp, fc !!soil moisture, wilting point and field capacity of current cell
REAL (KIND = float) :: actualE, actualT !!actual evaporation and transpiration of current cell
REAL (KIND = double) :: soilDepthCell !!soil depth of current cell [m]
REAL (KIND = float) :: prunoffcum !!cumulated runoff computed by Ross at previous timestep [cm]
REAL (KIND = float) :: pinfilcum !!cumulated infiltration computed by Ross at previous timestep [cm]
REAL (KIND = float) :: pdrncum !!cumulated drainage computed by Ross at previous timestep [cm]
REAL (KIND = float) :: pevapcum !!cumulated evapotranspiration used and possibly modified by Ross at previous timestep [cm]
REAL (KIND = float) :: water !!water in soil [m]
REAL (KIND = float) :: ddx	!actual cell length (m)
REAL (KIND = float) :: meanHydCond  !!mean hydraulic conductivity 
                             !for computing lateral subsurface flow (m/s)
REAL (KIND = float) :: ds !!thickness of the saturated depth (m)
REAL (KIND = float) :: cellWidth !!cell width (m)
REAL (KIND = float) :: ijSlope
LOGICAL             :: snowCovered, iceCovered


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

!reset variables for the time step 
QinSoilSub = 0.
QoutSoilSub = 0.
runoff = 0.
rainBalance = 0.


!update PET
!----------
CALL PETupdate (time)


!compute lateral subsurface intercell fluxes
!------------------------------------------
DO i = 1, mask % idim
    DO j = 1, mask % jdim
        IF ( mask % mat (i,j) == mask % nodata ) CYCLE
        
         !check balance id
         IF ( balanceId % mat (i,j) == LANDPLAIN .OR. &
              balanceId % mat (i,j) == LAKE ) THEN
             CYCLE
         END IF
         
        !set downstream cell  (is,js)
        CALL DownstreamCell(i, j, flowdir % mat(i,j), is, js, ddx, flowdir ) 
        
        IF ( CheckOutlet (i, j, is, js, flowdir) ) CYCLE
      
        !harmonic mean saturated conductivity
        meanHydCond = 2. * ksat_sub % mat (i,j) * ksat_sub % mat (is,js) / &
                    ( ksat_sub % mat (i,j) + ksat_sub % mat (is,js) )
      
        !thickness of the saturated depth
        ds = soilDepthTZ % mat(i,j) * soilMoistureTZ % mat (i,j) / thetas % mat (i,j)
        
        !cell width
        cellWidth = ( CellArea (soilMoisture, i, j) ) ** 0.5
        
        !local slope
        ijSlope = ( dem % mat (i,j) - dem % mat (is,js) ) / ddx
        
        IF ( ijSlope <= 0. ) THEN
            ijSlope = 0.0001
        END IF
      
        QoutSoilSub % mat (i,j) = cellWidth * ds * meanHydCond *  ijSlope
      
        !output becomes input of downstream cell
        QinSoilSub % mat (is,js) = QinSoilSub % mat (is,js) + QoutSoilSub % mat (i,j)
      
    END DO
END DO

!loop through mask
DO i = 1, mask % idim 
  DO j = 1, mask % jdim
     IF ( mask % mat (i,j) == mask % nodata ) CYCLE
        
     !reset fluxes of current cell
     runoffcell = 0.
     percolationcellRZ = 0.
     percolationcellTZ = 0.
     
     !set balance id
     id = balanceId % mat (i,j)
     
     !rain contributes to soil balance
     rainBalance % mat (i,j) = rain % mat (i,j)
     
     !add irrigation
     IF ( ALLOCATED ( irrigation % mat ) ) THEN
         rainBalance % mat (i,j) = rainBalance % mat (i,j) + irrigation % mat (i,j)
     END IF
    
     snowCovered = .FALSE.
     iceCovered  = .FALSE.
     !check if snow and ice cover the current cell
     IF ( dtIce > 0 ) THEN  !glaciers are simulated
        IF ( icewe % mat (i,j) <= 0. ) THEN !liquid water in ice 
                                            !contributes to soil balance
            rainBalance % mat (i,j) =  rainBalance % mat (i,j) + &
                                       waterInIce % mat (i,j) / dtSoilBalance
            waterInIce % mat (i,j) = 0.
            IF ( dtSnow > 0 ) THEN !snow is simulated
               IF ( swe % mat (i,j) <= 0. ) THEN !liquid water in snow
                                                 ! contributes to soil balance
                  rainBalance % mat (i,j) =  rainBalance % mat (i,j) + &
                                  waterInSnow % mat (i,j) / dtSoilBalance
                  waterInSnow % mat (i,j) = 0.
               ELSE
                  rainBalance % mat (i,j) = 0.
                  snowCovered = .TRUE.
               END IF
            END IF
        ELSE  !soil is frozen, no vertical fluxes in it 
            !runoff % mat(i,j) = 0.
            !infilt % mat(i,j) = 0.
            !percolation % mat(i,j) = 0.
            !et % mat(i,j) = 0.
            rainBalance % mat (i,j) = 0.
            iceCovered = .TRUE.
            
        END IF
     END IF
    
    IF ( dtSnow > 0 .AND. dtIce <= 0 ) THEN
       IF ( swe % mat (i,j) > 0.) THEN
          !soil is frozen, no vertical fluxes in it 
          !runoff % mat(i,j) = 0.
          !infilt % mat(i,j) = 0.
          !percolation % mat(i,j) = 0.
          !et % mat(i,j) = 0.
          rainBalance % mat (i,j) = 0.
          snowCovered = .TRUE.
       ELSE
          rainBalance % mat (i,j) = rainBalance % mat (i,j) + &
                                    waterInSnow % mat (i,j) / dtSoilBalance
          waterInSnow % mat (i,j) = 0.
       END IF
    END IF
      
   !decide if storm or interstorm period
    IF ( rainFlag % mat(i,j) == 0 .AND. rainBalance % mat (i, j) >= thresholdStartEvent ) THEN
			!storm starts
			rainFlag % mat(i,j) = 1
			interstormDuration % mat(i,j) = 0 !reset interstorm duration
            
			IF (infiltrationModel == SCS_CN) THEN
				CALL Sat2scn (curveNumber % mat(i,j), storativity % mat(i,j), &
					                soilSatRZ % mat(i,j), sEff % mat(i,j) )
			END IF
      
            IF (infiltrationModel == PHILIPEQ .OR. &
                  infiltrationModel == GREEN_AMPT) THEN
				        ism % mat(i,j) = soilMoistureRZ % mat(i,j)
            END IF
  
            IF (infiltrationModel == ROSS_BC .OR. &
                  infiltrationModel == ROSS_VG) THEN
				        runoffcum % mat(i,j) = 0.
                drncum % mat(i,j) = 0.
                infilcum % mat(i,j) = 0.
                evapcum % mat(i,j) = 0.
            END IF
     

    ELSE IF ( rainFlag % mat(i,j) == 1 .AND. rainBalance % mat (i, j) >= thresholdStartEvent ) THEN
			!storm continues
			rainFlag % mat(i,j) = 2
    ELSE IF ( rainBalance % mat (i, j) < thresholdStartEvent ) THEN
			!update interstorm duration
			interstormDuration % mat(i,j) = interstormDuration % mat(i,j) + dtSoilBalance
    END IF
	!check if storm event must be terminated
	IF (interstormDuration % mat(i,j) >= interstorm) THEN
          rainFlag % mat(i,j) = 0
            
          IF (infiltrationModel == SCS_CN) THEN
			      !reset cumulated rainfall
			      raincum % mat(i,j) = 0.
          END IF
          
          IF (infiltrationModel == PHILIPEQ .OR. &
              infiltrationModel == GREEN_AMPT ) THEN
			      !reset cumulated infiltration
			      cuminf % mat(i,j) = 0.
          END IF
    END IF
    
    !compute vertical fluxes according to soil balance computation scheme 
    SELECT CASE ( id )
    CASE (LAKE)
       !Lake drains subsurface flux
       !runoff can become < 0 when pet > rain + QinSoilSub
       runoffcell = rainBalance % mat (i, j) - pet % mat(i,j) + &
                  QinSoilSub % mat(i,j) / CellArea (QinSoilSub, i, j) 
       et % mat(i,j) = pet % mat(i,j) !et is potential 
       percolationcellRZ = 0. ! percolation is zero
       percolationcellTZ = 0.
       
       
    !CASE (LAKE_EWB)
    !    Tdeep = 25. ![°C]
        ! GR TODO
        !CALL ET_energy_nr_lago(SM%mat(iin,jin),albedo%mat(iin,jin),bpress%mat(iin,jin),&
		      !                              dtm%mat(iin,jin),temperatura_oss%mat(iin,jin),&
		      !                              radiazione_oss%mat(iin,jin),vento_oss%mat(iin,jin),umidita_oss%mat(iin,jin),&
		      !                              indice_nuv%mat(iin,jin),ET%mat(iin,jin),&
		      !                              Ts%mat(iin,jin),Rnetta%mat(iin,jin),Xle%mat(iin,jin),Hse%mat(iin,jin),&
		      !                              Ge%mat(iin,jin),Tdeep)                    	   
			!runoffcell = rainBalance % mat (i, j) - et % mat(i,j) + &
   !               QinSoilSub % mat(i,j) / CellArea (QinSoilSub, i, j)
   !        percolationcell = 0. !deep percolation is zero
           
    CASE DEFAULT ! all other cells(SLOPE, SLOPE_EWB, SLOPE_EWB_WET, CHANNEL, CHANNEL_EWB, &
         !CHANNEL_EWB_WET, PLAIN, PLAIN_EWB, PLAIN_EWB_WET)
      
         !compute infiltration and runoff during storm period
         IF ( rainFlag % mat(i,j) /= 0 ) THEN 
        
              SELECT CASE (infiltrationModel)  
          
              CASE (SCS_CN) !use soil conservation services curve number method
                infilt % mat (i,j) = SCScurveNumber (rain = rainBalance % mat (i, j), &
                                         raincum = raincum % mat (i,j), &
                                         c = abstractionRatio % mat(i,j), &
                                         s = sEff % mat(i,j), &
                                         dt = dtSoilBalance, runoff = runoffcell) 

          
              CASE (PHILIPEQ) !use Philips infiltration model  
                infilt % mat (i,j) = Philip (rain = rainBalance % mat (i, j), &
                                         psic = psic % mat(i,j), &
                                         ksat = ksat % mat(i,j), &
                                         theta = soilMoistureRZ % mat(i,j), &
                                         thetas = thetas % mat(i,j), &
                                         thetar = thetar % mat(i,j), &
                                         psdi = psdi % mat (i,j), &
                                         itheta = ism % mat (i,j), dt = dtSoilBalance, &
                                         cuminf = cuminf % mat (i,j) ) 
         
                runoffcell = rainBalance % mat (i, j) - infilt % mat (i,j)
         
          
              CASE (GREEN_AMPT) !use Green-Ampt infiltration model  
                infilt % mat (i,j) = GreenAmpt (rainFlag % mat(i,j), &
                                         rain = rainBalance % mat (i, j),  &
                                         ksat = ksat % mat(i,j), &
                                         theta = soilMoistureRZ % mat(i,j), &
                                         thetas = thetas % mat(i,j), &
                                         thetar = thetar % mat(i,j), &
                                         phy = phy % mat(i,j), &
                                         itheta = ism % mat (i,j), dt = dtSoilBalance, &
                                         cuminf = cuminf % mat (i,j) )    
         
                runoffcell = rainBalance % mat (i, j) - infilt % mat (i,j) 
         
             !CASE (ROSS_BC)
                 !the subroutine SolveRichardsBC is always called even during interstorm period  
         
              END SELECT
      ELSE !interstorm period, runoff is zero. Low intensity rainfall infiltrates
          runoffcell = 0.
          infilt % mat(i,j) = rainBalance % mat (i, j)
      END IF
      
      !compute actual evapotranspiration
      !---------------------------------
      !IF (id == SLOPE_EWB .OR. id == SLOPE_EWB_WET .OR. id == CHANNEL_EWB .OR. &
      !    id == CHANNEL_EWB_WET .OR. id == PLAIN_EWB .OR. id == PLAIN_EWB_WET) THEN
             
           !solve energy balance to compute actual evapotranspiration  
           !IF ( WetlandIsWet (wetland, iin, jin, time) ) THEN
       
               ! TODO
               !  Tsprof=26.5
				           !
		             !CALL ET_energy_nr_wetland(SM%mat(iin,jin),albedo%mat(iin,jin),bpress%mat(iin,jin),&
		             !                       dtm%mat(iin,jin),temperatura_oss%mat(iin,jin),&
		             !                       radiazione_oss%mat(iin,jin),vento_oss%mat(iin,jin),umidita_oss%mat(iin,jin),&
		             !                       indice_nuv%mat(iin,jin),ET%mat(iin,jin),&
		             !                       Ts%mat(iin,jin),Rnetta%mat(iin,jin),Xle%mat(iin,jin),Hse%mat(iin,jin),&
		             !                       Ge%mat(iin,jin),Tsprof)                    
		             !   
	              !  if(ET%mat(iin,jin)<0)then
		             !   ET%mat(iin,jin)=0
		             ! end if 
           !ELSE
         !      CALL ET_ENERGETICO_suolo(SM%mat(iin,jin),SMsat%mat(iin,jin),SMres%mat(iin,jin),albedo%mat(iin,jin),&
		       !                             res_stom_min%mat(iin,jin),LAI%mat(iin,jin),h_veg%mat(iin,jin),bpress%mat(iin,jin),&
		       !                             dtm%mat(iin,jin),B%mat(iin,jin),fraz_veg%mat(iin,jin),temperatura_oss%mat(iin,jin),&
		       !                             radiazione_oss%mat(iin,jin),vento_oss%mat(iin,jin),umidita_oss%mat(iin,jin),&
		       !                             indice_nuv%mat(iin,jin),WP%mat(iin,jin),FC%mat(iin,jin),ET%mat(iin,jin),&
		       !                             Ts%mat(iin,jin),Rnetta%mat(iin,jin),Xle%mat(iin,jin),Hse%mat(iin,jin),&
		       !                             Ge%mat(iin,jin),Ta_prec%mat(iin,jin),Ts_prec%mat(iin,jin))
		       !     Ta_prec%mat(iin,jin)=temperatura_oss%mat(iin,jin)
         !           Ts_prec%mat(iin,jin)=Ts%mat(iin,jin)
	        !        IF(ET%mat(iin,jin)<0)THEN
		       !        ET%mat(iin,jin)=0
		       !     END IF
		       !     IF (PRESENT(stato_neve)) THEN
					    !IF (stato_neve%mat(iin,jin) /= 0.0) THEN
					    !    ET % mat(iin,jin) = 0.0
					    !    IF(temperatura_oss%mat(iin,jin)<0.)THEN
		       !               Ts%mat(iin,jin)=temperatura_oss%mat(iin,jin)
	        !                ELSE
		       !               Ts%mat(iin,jin)=0.
	        !                END IF
         !                   Rnetta%mat(iin,jin) = 0.0
         !                   Xle%mat(iin,jin) = 0.0
		       !             Ge%mat(iin,jin) = 0.0
		       !             Hse%mat(iin,jin) = 0.0
   		    !            END IF    
           !END IF
      !ELSE 
              
           IF ( snowCovered .OR. iceCovered ) THEN
               et % mat(i,j) = 0.0
           ELSE
              !compute actual evaporation as a fraction alpha of potential
               sm = soilMoistureRZ % mat(i,j)
               wp = wiltingPoint % mat(i,j)
               fc = fieldCapacity % mat(i,j)
               alpha = 0.082 * sm + 9.173 * sm**2. - 9.815 * sm**3.
               IF (alpha > 1.) alpha = 1.
               actualE = pet % mat(i,j) * alpha
			    
               !compute actual transpiration as a fraction beta of potential		
               beta =  (sm - wp) / (fc - wp)
               IF (beta > 1.) beta = 1.
               IF (beta < 0.) beta = 0.
               actualT = pet % mat(i,j) * beta
				   
		      !compute actual evapotranspiration
		      et % mat(i,j) = vf % mat(i,j) * actualT + (1. - vf % mat(i,j)) * actualE
		      IF ( et % mat(i,j) < 0.) THEN
			      et % mat(i,j) = 0.0
		      END IF
         END IF
      !END IF
      
      !compute percolation and capilary rise
      !------------------------------------------------------------------------
      CALL PercolationAndCaprise (id, i, j, rainBalance % mat (i, j), vadose, pet, &
                         soilDepthCell, percolationcellRZ, &
                         percolationcellTZ, runoffcell)
   

      !update soil moisture
      !--------------------
      SELECT CASE (infiltrationModel)
      CASE (SCS_CN, PHILIPEQ, GREEN_AMPT)
          CALL UpdateSoilMoisture (id, i, j, dtSoilBalance, &
                                   soilDepthCell, &
                                   percolationcellRZ, &
                                   percolationcellTZ, &
                                   runoffcell,  &
                                   QinSoilSub % mat (i,j), &
                                   QoutSoilSub % mat (i,j) )   
          
         

      CASE (ROSS_BC)
          !store values of cumulated variables at previous time step
          pevapcum = evapcum % mat(i,j)
          prunoffcum = runoffcum % mat(i,j)
          pinfilcum = infilcum % mat(i,j)
          pdrncum = drncum % mat(i,j)
 
          !update soil moisture profile and cumulated fluxes in cm
          CALL SolveRichardsBC (idt = dtSoilBalance, &
                                rain = rainBalance % mat (i, j) , &
                                pet = et % mat (i,j), &
                                p = parameters (cell % mat(i,j) ), &
                                Scell = S (cell % mat(i,j),:), &
                                dxcell = dx (cell % mat(i,j),:), &
                                nsol = 0, &
                                h0 = h0grid % mat(i,j), &
                                nsteps = nsteps % mat(i,j) , &
                                evap = evapcum % mat(i,j), &
                                runoff = runoffcum % mat(i,j), &
                                infil = infilcum % mat(i,j), &
                                drn = drncum % mat(i,j), &
                                watercontent = water )
          !compute fluxes of this time step in m/s
          runoffcell = ( runoffcum % mat(i,j) - prunoffcum ) / 100. / dtSoilBalance
          IF (runoffcell < 0.) runoffcell = 0.

          infilt % mat (i,j)  = ( infilcum % mat(i,j) - pinfilcum ) / 100. / dtSoilBalance
          IF (infilt % mat (i,j)  < 0.) infilt % mat (i,j)  = 0.

          et % mat (i,j)  = ( evapcum % mat(i,j) - pevapcum ) / 100. / dtSoilBalance
          IF (et % mat (i,j)  < 0.) et % mat (i,j)  = 0.

          percolationcellTZ  = ( drncum % mat(i,j) - pdrncum ) / 100. / dtSoilBalance
          IF (percolationcellTZ  < 0.) percolationcellTZ  = 0.

          !compute average soil moisture
          soilMoisture % mat (i,j) = water / soilDepthCell

       END SELECT
              
    END SELECT

    !update percolation grid
    percolation % mat (i,j) = percolationcellTZ
    
    !update runoff grid
    runoff % mat (i,j) = runoffcell
     
  END DO
END DO
 
RETURN
END SUBROUTINE SolveSoilBalance
  
  
!==============================================================================
!| Description: 
!   return .TRUE. if wetland is flooded
FUNCTION WetlandIsWet   & 
  !
  (wetland,row,col,time)            &
  RESULT (wet)

IMPLICIT NONE
! function arguments 
! Scalar arguments with intent(in):

TYPE (grid_integer):: wetland
INTEGER, INTENT (IN) :: row, col
TYPE (DateTime), INTENT (IN) :: time

!Local scalar:
LOGICAL :: wet
INTEGER :: code
INTEGER :: i

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

IF (ALLOCATED (wetland % mat) ) THEN
    code = wetland % mat(row,col)
    wet = .FALSE.

    DO i = 1, SIZE(wetCode)
      IF ( code == wetCode (i) ) THEN
          IF (time >= wetBegin (i) .AND. time <= wetEnd (i) ) THEN
          wet = .TRUE.    
          EXIT
        END IF
      END IF
    END DO
ELSE
    wet = .FALSE.
END IF

END FUNCTION WetlandIsWet
  
  

  
!==============================================================================
!| Description: 
!   set initial condition for soilwater balance
SUBROUTINE SetInitialCondition  & 
!
  (iniDB)

IMPLICIT NONE

! Arguments with intent(in):

TYPE (IniList), INTENT(IN) :: iniDB

!Local declaration:
INTEGER (KIND = short) :: i, j
REAL (KIND = float) :: scalar

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

!mandatory variables

! root-zone soil saturation degree
IF (SectionIsPresent('saturation-rz', iniDB)) THEN
     IF (KeyIsPresent ('scalar', iniDB, 'saturation-rz') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'saturation-rz')
        CALL NewGrid (soilSatRZ, mask, scalar)
    ELSE
        CALL GridByIni (iniDB, soilSatRZ, section = 'saturation-rz')
    END IF
ELSE !grid is mandatory: stop the program if not present
   CALL Catch ('error', 'SoilBalance',   &
			   'error in loading saturation-rz: ' ,  &
			    argument = 'section not defined in ini file' )
END IF

! transmission-zone soil saturation degree
IF (SectionIsPresent('saturation-tz', iniDB)) THEN
     IF (KeyIsPresent ('scalar', iniDB, 'saturation-tz') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'saturation-tz')
        CALL NewGrid (soilSatTZ, mask, scalar)
    ELSE
        CALL GridByIni (iniDB, soilSatTZ, section = 'saturation-tz')
    END IF
ELSE !grid is mandatory: stop the program if not present
   CALL Catch ('error', 'SoilBalance',   &
			   'error in loading saturation-tz: ' ,  &
			    argument = 'section not defined in ini file' )
END IF

! allocate mean saturation map
CALL NewGrid (soilSat, mask, 0.)

!allocate and set soil moisture
CALL NewGrid (soilMoisture, mask)
CALL NewGrid (soilMoistureRZ, mask)
CALL NewGrid (soilMoistureTZ, mask)
        
DO j = 1, mask % jdim
    DO i = 1, mask % idim
        SELECT CASE ( balanceId % mat (i,j) )
            CASE(LAKE)  !lake cells are saturated by definition
				soilSat % mat(i,j) = 1.
                soilSatRZ % mat(i,j) = 1.
                soilSatTZ % mat(i,j) = 1.
				soilMoisture % mat(i,j) = 1.
                soilMoistureRZ % mat(i,j) = 1.
                soilMoistureTZ % mat(i,j) = 1.
			CASE DEFAULT 
				soilMoistureRZ % mat(i,j) = thetar % mat(i,j) + &
                                    soilSatRZ % mat(i,j) * &
                                    (thetas % mat(i,j) - &
                                    thetar % mat(i,j) )
                
                soilMoistureTZ % mat(i,j) = thetar % mat(i,j) + &
                                    soilSatTZ % mat(i,j) * &
                                    (thetas % mat(i,j) - &
                                    thetar % mat(i,j) )
                
				    
        END SELECT
    END DO
END DO



!optional variables:

! interstorm duration 
IF (SectionIsPresent('interstorm-duration', iniDB)) THEN
     IF (KeyIsPresent ('scalar', iniDB, 'interstorm-duration') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'interstorm-duration')
        CALL NewGrid (interstormDuration, mask, INT(scalar))
    ELSE
        CALL GridByIni (iniDB, interstormDuration, section = 'interstorm-duration')
    END IF
ELSE !grid is optional: set to default = 0
   CALL NewGrid ( interstormDuration, mask, 0 )
END IF


! precipitation status 
IF (SectionIsPresent('precipitation-status', iniDB)) THEN
     IF (KeyIsPresent ('scalar', iniDB, 'precipitation-status') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'precipitation-status')
        CALL NewGrid (rainFlag, mask, INT(scalar))
    ELSE
        CALL GridByIni (iniDB, rainFlag, section = 'precipitation-status')
    END IF
ELSE !grid is optional: set to default = 0
   CALL NewGrid ( rainFlag, mask, 0 )
END IF


! variables for SCS-CN method
IF ( infiltrationModel == SCS_CN ) THEN
    
   ! effective soil retention capacity of SCS-CN method [mm]
   IF (SectionIsPresent('soil-retention', iniDB)) THEN
        IF (KeyIsPresent ('scalar', iniDB, 'soil-retention') ) THEN
           scalar = IniReadReal ('scalar', iniDB, 'soil-retention')
           CALL NewGrid (sEff, mask, scalar )
        ELSE
           CALL GridByIni (iniDB, sEff, section = 'soil-retention')
        END IF
   ELSE !grid is optional: set to default = 0
       CALL NewGrid ( sEff, mask, 0. )
   END IF
   
   ! cumulative precipitation
   IF (SectionIsPresent('cumulative-precipitation', iniDB)) THEN
        IF (KeyIsPresent ('scalar', iniDB, 'cumulative-precipitation') ) THEN
           scalar = IniReadReal ('scalar', iniDB, 'cumulative-precipitation')
           CALL NewGrid (cuminf, mask, scalar )
        ELSE
           CALL GridByIni (iniDB, cuminf, section = 'cumulative-precipitation')
        END IF
   ELSE !grid is optional: set to default = 0
       CALL NewGrid ( cuminf, mask, 0. )
   END IF
   
END IF


! variables for Philip and Green-Ampt methods
IF ( infiltrationModel == PHILIPEQ .OR. &
     infiltrationModel == GREEN_AMPT ) THEN
    
   ! cumulative infiltration
   IF (SectionIsPresent('cumulative-infiltration', iniDB)) THEN
        IF (KeyIsPresent ('scalar', iniDB, 'cumulative-infiltration') ) THEN
           scalar = IniReadReal ('scalar', iniDB, 'cumulative-infiltration')
           CALL NewGrid (cuminf, mask, scalar )
        ELSE
           CALL GridByIni (iniDB, cuminf, section = 'cumulative-infiltration')
        END IF
   ELSE !grid is optional: set to default = 0
       CALL NewGrid ( cuminf, mask, 0. )
   END IF
END IF


!state variable initialization
!IF (SectionIsPresent('initial-saturation', iniDB)) THEN
!    
!    !cold start
!    IF ( KeyIsPresent ('cold', iniDB, section = 'initial-saturation') ) THEN
!        !allocate state variables
!        CALL NewGrid (sEff, mask, 0.)
!        CALL NewGrid (rainFlag, mask, 0)
!        CALL NewGrid (interstormDuration, mask, 0)
!          
!        !initial saturation
!        isd = IniReadReal ('cold', iniDB, section = 'initial-saturation')
!            
!        CALL Catch ('info', 'SoilBalance: ',     &
!                    'initial degree of saturation: ', &
!                    argument = ToString(isd))
!    	
!        !same value for root and transmission zones
!        CALL NewGrid (soilSat, mask, isd)
!        CALL NewGrid (soilSatRZ, mask, isd)
!        CALL NewGrid (soilSatTZ, mask, isd)
!          
!          
!        !allocate and set soil moisture
!        CALL NewGrid (soilMoisture, mask)
!        CALL NewGrid (soilMoistureRZ, mask)
!        CALL NewGrid (soilMoistureTZ, mask)
!        
!        DO i = 1, mask % idim
!            DO j = 1, mask % jdim
!                SELECT CASE ( balanceId % mat (i,j) )
!                    CASE(LAKE)
!					    soilSat % mat(i,j) = 1.
!                        soilSatRZ % mat(i,j) = 1.
!                        soilSatTZ % mat(i,j) = 1.
!					    soilMoisture % mat(i,j) = 1.
!				    CASE DEFAULT 
!					    soilMoisture % mat(i,j) = thetar % mat(i,j) + &
!                                            soilSat % mat(i,j) * &
!                                            (thetas % mat(i,j) - &
!                                            thetar % mat(i,j) )
!               
!                     !lake cells are saturated by definition
!				    
!                END SELECT
!            END DO
!        END DO
!        
!       !same initial soil mositure for root and transmission zones
!       soilMoistureRZ = soilMoisture
!       soilMoistureTZ = soilMoisture
!		              
!    ELSE     !hot start
!		!soil moisture
!	    !TODO HOT START FOR ROOT AND TRANSMISSION ZONES
!        CALL GridByIni (iniDB, soilMoisture, section = 'initial-saturation')
!        IF  ( .NOT. CRSisEqual (mask = mask, grid = soilMoisture, &
!                                checkCells = .TRUE.) ) THEN
!            CALL Catch ('error', 'SoilBalance ',   &
!		        'wrong spatial reference in soil-moisture' )
!        END IF
!            
!        !compute soil relative saturation
!        CALL NewGrid (soilSat, mask)
!        
!        DO i = 1, mask % idim
!            DO j = 1, mask % jdim
!                SELECT CASE ( balanceId % mat (i,j) )
!                    !lake cells are saturated by definition
!				    CASE(LAKE)
!					    soilSat % mat(i,j) = 1.
!					    soilMoisture % mat(i,j) = thetas % mat(i,j)
!                    CASE DEFAULT
!					    soilSat % mat(i,j) = ( soilMoisture % mat(i,j) - &
!                                            thetar % mat(i,j)) / & 
!                                          ( thetas % mat(i,j) - &
!                                            thetar % mat(i,j) )
!                END SELECT
!            END DO
!        END DO
!          
!    
!        ! effective soil retention capacity of SCS-CN method [mm]
!        IF (infiltrationModel == SCS_CN ) THEN
!		    IF (SectionIsPresent('soil-retention', iniDB)) THEN
!                CALL GridByIni (iniDB, sEff, section = 'soil-retention')
!                IF  ( .NOT. CRSisEqual (mask = mask, grid = sEff, &
!                                        checkCells = .TRUE.) ) THEN
!                  CALL Catch ('error', 'SoilBalance',   &
!			            'wrong spatial reference in soil retention capacity sEff' )
!                END IF
!            ELSE
!                CALL Catch ('error', 'SoilWaterBalance: ',   &
!			    'missing soil-retention section in configuration file' )
!            END IF
!        END IF
!        
!        
!        
!         !cumulative infiltration
!        IF (infiltrationModel == PHILIPEQ .OR. &
!            infiltrationModel == GREEN_AMPT ) THEN
!		    IF (SectionIsPresent('cumulative-infiltration', iniDB)) THEN
!                CALL GridByIni (iniDB, cuminf, section = 'cumulative-infiltration')
!                  IF  ( .NOT. CRSisEqual (mask = mask, grid = cuminf, &
!                                    checkCells = .TRUE.) ) THEN
!                   CALL Catch ('error', 'SoilBalance',   &
!			        'wrong spatial reference in cumulative infiltration' )
!                 END IF
!            ELSE
!                CALL Catch ('error', 'SoilWaterBalance: ',   &
!			    'missing cumulative-infiltration section in configuration file' )
!            END IF
!        END IF
!        
!            
!		!precipitation status
!		IF (SectionIsPresent('precipitation-status', iniDB)) THEN
!            CALL GridByIni (iniDB, rainFlag, section = 'precipitation-status')
!            IF  ( .NOT. CRSisEqual (mask = mask, grid = rainFlag, &
!                                    checkCells = .TRUE.) ) THEN
!              CALL Catch ('error', 'SoilBalance',   &
!			    'wrong spatial reference in precipitation status rainFlag' )
!            END IF
!        ELSE
!            CALL Catch ('error', 'SoilBalance: ',   &
!			    'missing precipitation-status section in configuration file' )
!        END IF
!            
!		
!            
!    !interstorm duration
!		IF (SectionIsPresent('interstorm-duration', iniDB)) THEN
!        CALL GridByIni (iniDB, interstormDuration, section = 'interstorm-duration')
!        IF  ( .NOT. CRSisEqual (mask = mask, grid = interstormDuration, &
!                                checkCells = .TRUE.) ) THEN
!          CALL Catch ('error', 'SoilBalance',   &
!			    'wrong spatial reference in interstorm duration' )
!        END IF
!    ELSE
!        CALL Catch ('error', 'SoilBalance: ',   &
!			'missing interstorm-duration section in configuration file' )
!    END IF
!    	
!	END IF !hot start
!ELSE
!  CALL Catch ('error', 'SoilBalance: ',   &
!			'missing initial-saturation section in configuration file' )
!END IF

RETURN
END SUBROUTINE SetInitialCondition

  
!==============================================================================
!| Description: 
!   set wetland for soilwater balance
SUBROUTINE SetWetland  & 
!
(iniDB, inifile)

IMPLICIT NONE

! Arguments with intent(in):

TYPE (IniList), INTENT(IN) :: iniDB
CHARACTER (LEN = *), INTENT(IN) :: inifile

!Local declaration:
TYPE (Table) :: wetlandTable
INTEGER (KIND = short) :: k

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

IF (SectionIsPresent('wetland', iniDB)) THEN
  
    CALL GridByIni (iniDB, wetland, section = 'wetland')

    !read table of flooded period  
    CALL TableNew ( inifile, wetlandTable )
        
    !allocate arrays
    ALLOCATE ( wetCode ( TableGetNrows(wetlandTable) ) )
    ALLOCATE ( wetBegin ( TableGetNrows(wetlandTable) ) )
    ALLOCATE ( wetEnd ( TableGetNrows(wetlandTable) ) )
        
    !fill in arrays 
    DO k = 1, TableGetNrows(wetlandTable) 
        wetCode (k) = StringToShort ( wetlandTable % col(1) % row(k) )
        wetBegin (k) = wetlandTable % col(2) % row(k) 
        wetEnd (k) = wetlandTable % col(3) % row(k) 
    END DO            
END IF

RETURN
END SUBROUTINE SetWetland

!==============================================================================
!| Description: 
!   compute percolation and capilalry rise
SUBROUTINE PercolationAndCaprise  & 
!
(id, i, j, rain, vadose, pet, soilDepthCell, percolationcellRZ, &
    percolationcellTZ, runoffcell)

IMPLICIT NONE

! Arguments with intent(in):
INTEGER (KIND = short), INTENT(IN) :: id !!soil balance id
INTEGER (KIND = short), INTENT(IN) :: i,j
REAL (KIND = float), INTENT(IN) :: rain
TYPE (grid_real), INTENT(IN) :: vadose !!vadose zone depth
TYPE (grid_real), INTENT(IN) :: pet !!potential evapotranspiration

!Arguments with intent(out):
REAL (KIND = double), INTENT(OUT) :: soilDepthCell
REAL (KIND = float), INTENT(OUT) :: percolationcellRZ
REAL (KIND = float), INTENT(OUT) :: percolationcellTZ
REAL (KIND = double), INTENT(OUT) :: runoffcell

!Local declaration:
REAL (KIND = float) :: smAdjustedRZ !! used to prevent soil moisture going to zero
REAL (KIND = float) :: smAdjustedTZ !! used to prevent soil moisture going to zero
REAL (KIND = float) :: conductivityCellRZ !!partial saturation soil hydraulic 
                                        !conductivity of current cell [m/s]
REAL (KIND = float) :: conductivityCellTZ !!partial saturation soil hydraulic 
                                        !conductivity of current cell [m/s]
                                        !conductivity of current cell [m/s]
REAL (KIND = float) :: psiCell !!matric potential of current cell [m]
REAL (KIND = float) :: meanHydCond !!mean hydraulic conductivity used to 
                                   !compute capillary rise [m/s]
REAL (KIND = float) :: dsdt
!------------end of declaration------------------------------------------------

!calculate unsaturated hydraulic conductivity [m/s]

!root zone
IF (soilMoistureRZ % mat(i,j) <= thetar % mat(i,j)) THEN
    smAdjustedRZ = thetar % mat(i,j) + 0.001
ELSE
    smAdjustedRZ = soilMoistureRZ % mat(i,j) 
END IF

conductivityCellRZ = UnsHydCond (ksat = ksat % mat(i,j), &
                    theta = smAdjustedRZ, &
                    thetas = thetas % mat(i,j), &
                    thetar = thetar % mat(i,j), &
                    psdi = psdi % mat(i,j) ) 

!transmission zone
IF (soilMoistureTZ % mat(i,j) <= thetar % mat(i,j)) THEN
    smAdjustedTZ = thetar % mat(i,j) + 0.001
ELSE
    smAdjustedTZ = soilMoistureTZ % mat(i,j) 
END IF

conductivityCellTZ = UnsHydCond (ksat = ksat % mat(i,j), &
                    theta = smAdjustedTZ, &
                    thetas = thetas % mat(i,j), &
                    thetar = thetar % mat(i,j), &
                    psdi = psdi % mat(i,j) ) 

IF (id == LANDPLAIN ) THEN

         
    !set local soil depth
    !soilDepthCell = MIN (soildepth % mat(i,j), vadose % mat(i,j) )
    soilDepthCell = soildepth % mat(i,j)
    
    !compute root zone soil suction (m)
    psiCell = Psi (psic = psic % mat(i,j), theta = smAdjustedRZ, &
                    thetas = thetas % mat(i,j), &
                    thetar = thetar % mat(i,j), &
                    psdi = psdi % mat(i,j) ) 
				   
    !Compute harmonic mean among saturated conductivity 
    !at the interface with groundwater table and unsaturated 
    !conductivity of vadose zone
    meanHydCond = 2. * conductivityCellRZ * ksat % mat(i,j) / &
                  (conductivityCellRZ + Ksat % mat(i,j))
        

    !compute capillary rise and percolation 
    IF  ( vadose % mat (i,j) < soilDepthRZ % mat (i,j) )  THEN
        ! water table depth lies within the root zone or above ground
        saturatedByGroundwater = .TRUE.
    ELSE
        capRise % mat (i,j) =  meanHydCond * psiCell / vadose % mat (i,j)
        dsdt = ( soilMoistureRZ % mat(i,j) - thetar % mat(i,j) ) * &
                      soilDepthRZ % mat(i,j) / dtSoilBalance
        percolationcellRZ = MIN ( conductivityCellRZ, dsdt)
            
        dsdt = ( soilMoistureTZ % mat(i,j) - thetar % mat(i,j) ) * &
                      soilDepthTZ % mat(i,j) / dtSoilBalance
        percolationcellTZ = MIN ( conductivityCellTZ, dsdt) * &
                                percolationFactor % mat (i,j)
        
        saturatedByGroundwater = .FALSE.
    END IF
       
				    
ELSE !slope or channel cell
			soilDepthCell = soilDepth % mat(i,j)
			capRise % mat(i,j) = 0.
            dsdt = ( soilMoistureRZ % mat(i,j) - thetar % mat(i,j) ) * &
                      soilDepthRZ % mat(i,j) / dtSoilBalance
            percolationcellRZ = MIN ( conductivityCellRZ, dsdt)
            
            dsdt = ( soilMoistureTZ % mat(i,j) - thetar % mat(i,j) ) * &
                      soilDepthTZ % mat(i,j) / dtSoilBalance
            percolationcellTZ = MIN ( ksat % mat(i,j), dsdt) * &
                                percolationFactor % mat (i,j)
END IF

 
RETURN
END SUBROUTINE PercolationAndCaprise


!==============================================================================
!| Description: 
!   update soil moisture
SUBROUTINE UpdateSoilMoisture  & 
!
(id, i, j, dt, soilDepthCell, percolationcellRZ, percolationcellTZ, &
 runoffcell, Qin, Qout )

IMPLICIT NONE

! Arguments with intent(in):
INTEGER (KIND = short), INTENT(IN) :: id !!soil balance id
INTEGER (KIND = short), INTENT(IN) :: i,j
INTEGER (KIND = short), INTENT(IN) :: dt
REAL (KIND = double), INTENT(IN) :: soilDepthCell
REAL (KIND = float), INTENT(IN) :: Qin
REAL (KIND = float), INTENT(IN) :: Qout



!Arguments with intent(inout):
REAL (KIND = double), INTENT(INOUT) :: runoffcell
REAL (KIND = float), INTENT(INOUT) :: percolationcellRZ
REAL (KIND = float), INTENT(INOUT) :: percolationcellTZ

!Local declaration:
REAL (KIND = double) :: smPrevRZ !! root zone soil moisture of the previous time step
REAL (KIND = double) :: smPrevTZ !! transmission zone soil moisture of the previous time step
REAL (KIND = float)  :: check
REAL (KIND = double) :: dsRZ !!root zone soil mositure excess rate (m/s)
REAL (KIND = double) :: dsTZ !!transmission zone soil mositure excess rate (m/s)
REAL (KIND = float)  :: sdRZ  !!root zone soil depth (m)
REAL (KIND = float)  :: sdTZ  !!transmission zone soil depth (m)
REAL (KIND = double) :: subIn, subOut  !!subsurface fluxes
REAL (KIND = double) :: area !!cell area (m2)
!------------end of declaration------------------------------------------------

IF ( id == LANDPLAIN .AND. saturatedByGroundwater ) THEN
    ! water table depth lies within the root zone or above ground
    percolationcellRZ = 0.
    percolationcellTZ = 0.
    infilt % mat(i,j) = 0.
    et % mat(i,j) = pet % mat(i,j)
    capRise % mat (i,j) = et % mat(i,j)
    runoffcell = rainBalance % mat (i, j)
    soilMoistureRZ % mat(i,j) = thetas % mat (i,j)
    soilMoistureTZ % mat(i,j) = thetas % mat (i,j)
    soilMoisture % mat(i,j) = thetas % mat (i,j)
    balanceError % mat (i,j) = 0.
    RETURN
END IF

IF ( id == LAKE ) THEN
    soilMoistureRZ % mat(i,j) = thetas % mat (i,j)
    soilMoistureTZ % mat(i,j) = thetas % mat (i,j)
    soilMoisture % mat(i,j) = thetas % mat (i,j)
    balanceError % mat (i,j) = 0.
    deltaSoilMoisture % mat (i,j) = 0.
    RETURN
END IF

smPrevRZ = soilMoistureRZ % mat(i,j)
smPrevTZ = soilMoistureTZ % mat(i,j)

sdRZ = soilDepthRZ % mat (i,j)
sdTZ = soilDepthTZ % mat (i,j)

!IF ( soilDepthCell < soilDepthRZ % mat(i,j) ) THEN
!    sdRZ = soilDepthCell
!    sdTZ = 0.
!ELSE
!    sdRZ = soilDepthRZ % mat(i,j)
!    sdTZ = MAX (0.1, ( soilDepthCell - sdRZ ) ) 
!END IF

IF ( id == HILLSLOPE .OR. id == CHANNEL ) THEN
    area = CellArea (soilMoisture, i, j) 
    subIn = Qin / area
    subOut = Qout / area
ELSE IF ( id == LANDPLAIN ) THEN
    subIn = 0.
    subOut = 0.
END IF


IF (sdTZ > 0.) THEN
	soilMoistureTZ % mat(i,j) = smPrevTZ + & 
                        ( percolationcellRZ -  &
                        percolationcellTZ + &
                        subIn - subOut      )  /  &
                        sdTZ * dt
        
    IF (soilMoistureTZ % mat(i,j) < thetar % mat(i,j)) THEN
        soilMoistureTZ % mat(i,j) = thetar % mat(i,j)
    END IF
        
    IF ( soilMoistureTZ % mat(i,j) >  thetas % mat(i,j) ) THEN !saturation excess
        dsTZ = ( soilMoistureTZ % mat(i,j) - thetas % mat(i,j) ) * sdTZ / dt
            soilMoistureTZ % mat(i,j) = thetas % mat(i,j)
    ELSE
        dsTZ = 0.
    END IF
              
END IF

IF (sdRZ > 0.) THEN
	soilMoistureRZ % mat(i,j) = smPrevRZ + ( infilt % mat(i,j) - percolationcellRZ - &
                            et % mat(i,j) + capRise % mat(i,j) + dsTZ  )  / sdRZ * dt		   
              
END IF


!particular cases
!case 1: soil is too dry, soil moisture goes below residual value
IF (soilMoistureRZ % mat(i,j) < thetar % mat(i,j) ) THEN
  !compute soil moisture deficit (negative) converted to flux (m/s)
  dsRZ = ( soilMoistureRZ % mat(i,j) - thetar % mat(i,j) ) * soilDepthCell / dt
  
  !adjust evapotranspiration
  et % mat(i,j) = et % mat(i,j) + dsRZ
     
  !check for negative et
  IF (et % mat(i,j) < 0. ) et % mat(i,j) = 0.
     
  !limit soil moisture
  soilMoistureRZ % mat(i,j)  = thetar % mat(i,j)
  
  
!case 2: soil is too wet, soil moisture goes above saturated value  
ELSE IF ( soilMoistureRZ % mat(i,j) > thetas % mat(i,j) ) THEN
  !compute soil moisture excess converted to flux (m/s)
	dsRZ = ( soilMoistureRZ % mat(i,j) - thetas % mat(i,j) ) * sdRZ / dt
	!adjust runoff for saturation excess runoff and infiltration rate
	runoffcell = runoffcell + dsRZ
	infilt % mat(i,j) = infilt % mat(i,j) - dsRZ
    
    !limit soil moisture to saturation value
    soilMoistureRZ % mat(i,j)  = thetas % mat(i,j)
  
    !check if infiltration is negative after correction
    !adjust capillary rise if present (plain cell)
	IF (infilt % mat(i,j) < 0.) THEN
      dsRZ = - infilt % mat(i,j)
      infilt % mat(i,j) = 0.
      IF (id == LANDPLAIN ) THEN
         capRise % mat(i,j) = capRise % mat(i,j) - dsRZ
      END IF
	END IF
	

END IF

!balance error
check = infilt % mat(i,j) +   &
        capRise % mat(i,j) -  &
        percolationcellTZ  -    &
        et % mat(i,j) +       &
        subIn - subOut - &
        (soilMoistureRZ % mat(i,j) - smPrevRZ) * sdRZ / dt - &
        (soilMoistureTZ % mat(i,j) - smPrevTZ) * sdTZ / dt

balanceError % mat (i,j) = check * dt * 1000.


!soil mositure variation
deltaSoilMoisture % mat (i,j) = (soilMoistureRZ % mat(i,j) - smPrevRZ) * sdRZ + &
        (soilMoistureTZ % mat(i,j) - smPrevTZ) * sdTZ 


!update relative saturation and mean soil moisture
soilSatRZ % mat(i,j) = ( soilMoistureRZ % mat(i,j) - thetar % mat(i,j) ) / &
                     ( thetas % mat(i,j) - thetar % mat(i,j) )
soilSatTZ % mat(i,j) = ( soilMoistureTZ % mat(i,j) - thetar % mat(i,j) ) / &
                     ( thetas % mat(i,j) - thetar % mat(i,j) )
soilSat % mat(i,j) = ( soilSatRZ % mat(i,j) * sdRZ  + &
                       soilSatTZ % mat(i,j) * sdTZ  ) / (sdRZ + sdTZ)
IF ( sdRZ + sdTZ > 0. ) THEN
     soilMoisture % mat(i,j) = ( soilMoistureRZ % mat(i,j) * sdRZ + &
                            soilMoistureTZ % mat(i,j) * sdTZ ) / (sdRZ + sdTZ)
ELSE
    soilMoisture % mat(i,j) = thetas % mat (i,j) 
END IF

RETURN
END SUBROUTINE UpdateSoilMoisture


END MODULE SoilBalance