!! Perform hydrological simulation
!|author:  <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a>
! license: <a href="http://www.gnu.org/licenses/">GPL</a>
!
!```
!  ________  _______       ________       _________    
! |\  _____\|\  ___ \     |\   ____\     |\___   ___\  
! \ \  \__/ \ \   __/|    \ \  \___|_    \|___ \  \_|  
!  \ \   __\ \ \  \_|/__   \ \_____  \        \ \  \   
!   \ \  \_|  \ \  \_|\ \   \|____|\  \        \ \  \  
!    \ \__\    \ \_______\    ____\_\  \        \ \__\ 
!     \|__|     \|_______|   |\_________\        \|__| 
!                            \|_________|           
!          S P A T I A L L Y - D I S T R I B U T E D 
!              H Y D R O L O G I C A L - M O D E L    
!```

!### Description 
! _FeST_ is a spatially distributed physically based hydrological
! model. _FeST_ is the acronym of “flash–Flood Event–based Spatially 
! distributed rainfall–runoff Transformation”. It was initially developed 
! by Mancini (1990), as a model oriented to the simulation of 
! rainfall-runoff transformation of single flood events. 
! Later the _FeST_ model was merged with the soil water balance scheme
! from TOPLATS model (Famiglietti and Wood, 1994), transforming it into a 
! continuous model (Montaldo et al., 2007). Then the _FeST_ code was 
! redesigned and rewritten from scratch while keeping the basic assumptions
! of the previous release (Rabuffetti et al., 2008). In 2011 the _FeST_
! was upgraded with a routine to solve the system of water mass and 
! energy balance in order to better simulate the actual evapotranspiration
! and interface the model to remotely sensed data (Corbari et al., 2011;
! Corbari & Mancini, 2014). At the same year, 2011, a new module for 
! simulating groundwater flux and river-groundwater interaction
! was developed and implemented in the _FeST_ (Ravazzani et al., 2011).
! In 2013 a new version of the code was released built on top of
! the MOSAICO library (Ravazzani, 2013). In 2014 the _FeST_ 
! was upgraded with a module for glaciers modelling (Boscarello et al., 2014).
! In 2021 a forest growth component was implemented in the _FeST_ 
! (Feki et al., 2021).
!
!
!### References
!
! Boscarello, L., Ravazzani, G., Rabuffetti, D., & Mancini, M.. (2014). 
!    Integrating glaciers raster-based modelling in large catchments 
!    hydrological balance: the Rhone case study. 
!    Hydrological processes, 28, 496–508.
!
!
! Corbari, C., & Mancini, M. (2014). Calibration and validation of a distributed 
!    energy–water balance model using satellite data of land surface temperature 
!    and ground discharge measurements. Journal of hydrometeorology, 15(1), 376-392.
!
! Corbari, C., Ravazzani, G., & Mancini, M. (2011). A distributed thermodynamic 
!   model for energy and mass balance computation: FEST–EWB. 
!   Hydrological Processes, 25(9), 1443-1452.
!
! Famiglietti, J.S., & Wood, E.F. (1994). Multiscale modeling of spatially 
!    variable water and energy balanceprocess. Water Resour. Res., 
!    30, 3061–3078, https://doi.org/10.1029/94WR01498. 
!
! Feki, M., Ravazzani, G., Ceppi, A., Pellicone, G., & Caloiero, T.. (2021). 
!  Integration of forest growth component in the fest-wb distributed 
!  hydrological model: the bonis catchment case study. Forests, 12(12).
!
! Mancini, M.. (1990). La modellazione distribuita della risposta idrologica: 
!     effetti della variabilità spaziale e della scala di rappresentazione 
!     del fenomeno dell'assorbimento. PhD thesis. Politecnico di Milano, 
!     Istituto di idraulica
!
! Montaldo, N., Ravazzani, G., & Mancini, M.. (2007). On the 
!    prediction of the toce alpine basin floods with distributed 
!    hydrologic models. Hydrological processes, 21, 608–621
!
! Rabuffetti, D., Ravazzani, G., Corbari, C., & Mancini, M.. (2008). 
!   Verification of operational quantitative discharge forecast (QDF) 
!   for a regional warning system – the AMPHORE case studies in the 
!   upper Po river. Natural hazards and earth system sciences, 8, 161–173.
!
! Ravazzani, G.. (2013). Mosaico, a library for raster based hydrological 
!   applications. Computers & geosciences, 51, 1–6.
!
! Ravazzani, G., Rametta, D., & Mancini, M.. (2011). Macroscopic 
!   cellular automata for groundwater modelling: a first approach. 
!   Environmental modelling & software, 26(5), 634–643.
!
!
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
PROGRAM Fest

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

USE HydroNetwork, ONLY : &
  !imported definitions:
  !NetworkSurfaceFlow, 
  NetworkGroundwater

USE Reservoirs, ONLY : &
  !imported definitions:
  Reservoir, &
  !Imported variables:
  dtReservoir, &
  !Imported routines:
  ReservoirSaveStatus


USE Diversions, ONLY : &
    !imported variables:
    dtDiversion, &
    !imported routines:
    DiversionSaveStatus

USE SpatialAverage, ONLY : &
   !imported routines:
   InitSpatialAverageMeteo, &
   InitSpatialAverageBalance, &
   InitSpatialAverageSnow, &
   InitSpatialAverageGlaciers, &
   InitSpatialAverageSediment, &
   InitSpatialAverageCanopy, &
   InitSpatialAveragePlants, &
   ComputeSpatialAverageMeteo,&
   ComputeSpatialAverageBalance, &
   ComputeSpatialAverageSnow,  &
   ComputeSpatialAverageGlaciers, &
   ComputeSpatialAverageSediment, &
   ComputeSpatialAverageCanopy, &
   ComputeSpatialAveragePlants, &
   ExportSpatialAverageMeteo, &
   ExportSpatialAverageBalance, &
   ExportSpatialAverageSnow, &
   ExportSpatialAverageGlaciers, &
   ExportSpatialAverageSediment, &
   ExportSpatialAverageCanopy, &
   ExportSpatialAveragePlants, &
   timeSpatialAverageMeteo, &
   timeSpatialAverageBalance,&
   timeSpatialAverageSnow, &
   timeSpatialAverageIce, &
   timeSpatialAverageSediment, &
   timeSpatialAverageCanopy, &
   timeSpatialAveragePlants

USE BasinAverage, ONLY : &
  !Imported routines:
  InitBasinAverage, &
  ReadPointFileBasinAverage, &
  ExportBasinAverage, &
  !Imported variables:
  dtBasinAverage

USE RasterExport, ONLY : &
    !Imported routines:
    InitRasterExport, &
    ExportMaps

USE Snow, ONLY : &
  !imported routines:
  SnowInit, SnowUpdate, &
  SnowPointInit, &
  !Imported variables:
  dtSnow, meltCoefficient, &
  swe, waterInSnow, snowMelt, &
  rainfallRate

USE Glacier, ONLY : &
  !imported routines:
  GlacierInit, GlacierPointInit, &
  GlacierUpdate, &
  !imported variables:
  dtIce, icewe, waterInIce, iceMelt

USE SoilBalance, ONLY : &
  !imported routines
  InitSoilBalance, SolveSoilBalance, &
  !imported variables:
  soilMoisture, soilMoistureRZ, &
  soilMoistureTZ, soilSatRZ, &
  soilSatTZ, sEff, rainFlag, &
  runoff, infilt, percolation, caprise, &
  soilDepth, wiltingPoint, fieldCapacity, &
  psdi, ksat, et, pet, dtSoilBalance, &
  QinSoilSub, soilSat, balanceError, &
  interstormDuration, deltaSoilMoisture

USE Infiltration, ONLY : &
  !imported variables:
  infiltrationModel, SCS_CN, &
  PHILIPEQ, GREEN_AMPT, &
  sEff, raincum, cuminf

USE Groundwater, ONLY : &
  !Imported routines:
  GroundwaterInit, &
  GroundwaterUpdate, &
  GroundwaterPointInit, &
  GroundwaterOutputInit, &
  GroundwaterOutput, &
  GroundwaterRiverUpdate, &
  !imported variables:
  dtGroundwater, &
  vadoseDepth, &
  riverGroundwaterInteract, &
  riverToGroundwater, &
  groundwaterToRiver

USE DischargeRouting, ONLY : &
  !imported routines:
  InitDischargeRouting, DischargeRoute, &
  DischargePointInit, &
  !imported variables:
  Qin, Qout, Pin, Pout, Plat, &
  dtDischargeRouting, &
  topWidth, waterDepth

USE Irrigation, ONLY : &
    !Imported routines:
    IrrigationConfig, IrrigationUpdate, &
    !Imported variables:
    dtIrrigation, &
    irrigationFlux

USE Sediment, ONLY: &
    !Imported routines:
    SedimentInit, InterrillDetachmentRate, &
    SedimentRouting, &
    !imported variables:
    interrillErosion, routeSediment, QoutSed, deltaSed

USE DomainProperties, ONLY: &
    !imported variables:
    mask, albedoGround, albedo, &
    albedo_loaded,  &
    !imported routines:
    DomainInit

USE MorphologicalProperties, ONLY: &
    !imported variables:
    dem, dem_loaded, &
    flowDirection, &
    streamNetwork, &
    !imported routines:
    MorphologyInit


USE Meteo, ONLY: &
  !Imported routines:
  MeteoInit, MeteoRead , &
  MeteoPointInit, &
  !Imported variables:
  dtMeteo

USE Precipitation, ONLY: &
   !imported variables:
   precipitationRate

USE PrecipitationDaily, ONLY: &
    !imported variables:
    precipitationRateDaily

USE AirTemperature, ONLY: &
    !imported variables:
    temperatureAir

USE AirTemperatureDailyMean, ONLY: &
    !imported variables:
    temperatureAirDailyMean

USE AirTemperatureDailyMax, ONLY: &
    !imported variables:
    temperatureAirDailyMax

USE AirTemperatureDailyMin, ONLY: &
    !imported variables:
    temperatureAirDailyMin

USE SolarRadiation, ONLY: &
    !imported variables:
    radiation, netRadiation
 
USE AirRelativeHumidity, ONLY: &
    !imported variables:
    relHumidityAir

USE WindFlux, ONLY: &
    !imported variables:
    windSpeed


USE Plants, ONLY : &
    !imported routines:
    PlantsConfig, PlantsGrow, &
    PlantsParameterUpdate, &
    !imported variables:
    dtPlants, lai, fvcover, &
    gpp, npp, carbonroot, &
    carbonstem, carbonleaf, dbh, &
    plantsHeight, density, stemyield, &
    rsMinLoaded, updatePlantsParameters

USE PlantsInterception, ONLY: &
    !imported variables:
    dtCanopyInterception, canopyStorage, &
    canopyPT, &
    !imported routines:
    InterceptionInit, Throughfall, &
    AdjustPT


USE DataTypeSizes, ONLY: &
    !Imported parameters:
    short, long, float

USE LogLib, ONLY: &
    !Imported routines:
    LogInit, LogStop, Catch, &
    !Imported variables:
    verbose

USE GeoLib, ONLY: &
    !Imported routines:
    GeoInit

USE Chronos, ONLY: &
    !Imported routines:
    GetMinute, GetHour, GetMonth, &
    GetDay, GetDayOfWeek, &
    !Imported definitions:
    DateTime, &
    !Imported variables:
    timeString, &
    !Imported operators:
    ASSIGNMENT (=), &
    OPERATOR (+), &
    OPERATOR (-), &
    OPERATOR (==), &
    OPERATOR (>), &
    OPERATOR (<), &
    OPERATOR (<=)

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

USE GridLib, ONLY : &
    !imported definition:
    grid_real, &
    !imported routines:
    NewGrid, ExportGrid, &
    !imported parameters:
    ESRI_BINARY

USE StringManipulation, ONLY : &
    !imported routines
    StringCompact, StringTokenize, &
    StringToLong

USE CronSchedule, ONLY : &
    !Imported types:
    CronTab, &
    !Imported routines:
    CronParseString, &
    CronIsTime

IMPLICIT NONE

!configuration
TYPE (IniList) :: festIni
CHARACTER (len = 300)  :: arg !! command line arguments
CHARACTER (len = 1000) :: string
LOGICAL  :: iniFound !!says if configuration file has been passed as command argument

!simulation time
TYPE (DateTime) :: time  !current step time
TYPE (DateTime) :: timeStart !start time 
TYPE (DateTime) :: timeStop !stop time

!meteo
TYPE (DateTime)        :: timeNewMeteo
TYPE (DateTime)        :: timeOutMeteo
INTEGER (KIND = short) :: dtOutMeteo 

!irrigation
TYPE (DateTime)        :: timeNewIrrigation
TYPE (DateTime)        :: timeOutIrrigation
INTEGER (KIND = short) :: dtOutIrrigation

!snow
TYPE (DateTime) :: timeNewSnow
TYPE (DateTime) :: timeOutSnow
INTEGER (KIND = short) :: dtOutSnow

!glacier
TYPE(DateTime) :: timeNewIce
TYPE(DateTime) :: timeOutIce
INTEGER (KIND = short) :: dtOutIce

!canopy interception
TYPE(DateTime) :: timeNewCanopyInterception
TYPE(grid_real) :: raineff ! throughfall[m/s]
INTEGER (KIND = short) :: dtOutCanopy
TYPE (DateTime) :: timeOutCanopy

!forest dynamic
TYPE(DateTime) :: timeNewPlants
INTEGER (KIND = short) :: dtOutPlants
TYPE (DateTime) :: timeOutPlants    
    
    
!soil balance
TYPE(DateTime) timeNewSoilBalance

!bilancio energetico
TYPE (grid_real)Ts !temperatura supeficiale
TYPE (grid_real)Rnetta
TYPE (grid_real)Xle
TYPE (grid_real)Hse
TYPE (grid_real)Ge
TYPE (grid_real)Ta_prec       !temperatura aria istante precedente
TYPE (grid_real)Ts_prec     !temperatura suolo istante precedente
	
!output
TYPE(DateTime) :: timeOutSoilBalance
INTEGER :: dtOutSoilBalance

! groundwater
INTEGER (KIND = short) :: dtOutGroundwater
TYPE(DateTime)         :: timeNewGroundwater
TYPE(DateTime)         :: timeOutGroundwater

	
!sediment
INTEGER               :: dtCalcSediment
TYPE(DateTime)        :: timeUpdateSediment
INTEGER               :: dtOutSediment
TYPE(DateTime)        :: timeOutSediment
!sediment routing
INTEGER               :: dtCalcSedimentRouting
TYPE(DateTime)        :: timeUpdateSedimentRouting
INTEGER               :: dtOutSedimentRouting
TYPE(DateTime)        :: timeOutSedimentRouting
!	TYPE(rete_osservativa)  :: xsOutSedimentRouting

! discharge routing
TYPE (DateTime) :: timeNewDischargeRouting
TYPE(grid_real)Velocita
TYPE (grid_real) storageChannelSurf

!basin average
TYPE (DateTime)        :: timeNewBasinAverage
INTEGER (KIND = short) :: dtOutBasinAverage

!raster export
TYPE (DateTime)        :: timeNewRasterExport

!hotstart
CHARACTER (len=1000)   :: pathHotstart
LOGICAL                :: saveLast
LOGICAL                :: hotstart
CHARACTER (LEN = 300)  :: cron
TYPE (CronTab)         :: cronHotstart

!general variables
INTEGER (KIND = short)  :: dtCalc  !calculation time step
INTEGER (KIND = short)  :: dtArray (4)
CHARACTER (len = 1000)  :: path_out

INTEGER :: k,i,j
LOGICAL :: fatta_previ = .FALSE.


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

!   splash screen
IF ( verbose ) THEN
   CALL Splash ()
END IF
    
!------------------------------------------------------------------------------
!                               initialize log
!------------------------------------------------------------------------------
 
CALL LogInit

!------------------------------------------------------------------------------
!                            initialize cartographic engine
!------------------------------------------------------------------------------

CALL GeoInit ('GeoLib.ini')

!------------------------------------------------------------------------------
!           check command line arguments to read configuration file name
!------------------------------------------------------------------------------

i = 1
iniFound = .FALSE.
DO WHILE ( .not. (arg == '') )

	CALL Getarg ( i, arg )
	SELECT CASE (arg)
	  CASE ( '-inifile' )
        i = i + 1
		CALL Getarg ( i , string )
		iniFound = .TRUE.
      CASE DEFAULT
			!case unknown
	END SELECT
	i = i + 1		
END DO

!------------------------------------------------------------------------------
!                 read name of configuration file if not yet defined
!------------------------------------------------------------------------------

IF ( .NOT. iniFound) THEN

  WRITE(*,*) 'Insert the name of configuration file and press enter'
  READ(*,*) string

END IF

CALL Catch ('info', 'Fest', &
         'start configuration from file: ', argument = string )

!------------------------------------------------------------------------------
!                                read configuration file
!------------------------------------------------------------------------------
IF ( verbose ) THEN
    WRITE(*,*)
    WRITE(*,*) 'Loading configuration...'
END IF

string = "./conf-test/fest.ini"

CALL IniOpen (string, festIni)

!------------------------------------------------------------------------------
!                               configure time
!------------------------------------------------------------------------------
CALL Catch ('info', 'Fest', 'time configuration' )

IF ( verbose ) THEN
    WRITE(*,*)
    WRITE(*,*) '...configuring time...'
END IF

!check whether time is present
IF ( .NOT. SectionIsPresent ('time', festIni) ) THEN  !section is mandatory
    CALL Catch ('error', 'Fest', 'time section missing in configuration file' )
END IF

!start time
IF (KeyIsPresent('start', festIni, section = 'time')) THEN	
   timeString = IniReadString ('start', festIni, section = 'time')
ELSE
    CALL Catch ('error', 'Fest', 'start missing in time' )
END IF

timeStart = timeString

CALL Catch ('info', 'Fest', 'start time: ', argument = timeString )

!stop time
IF (KeyIsPresent('stop', festIni, section = 'time')) THEN	
   timeString = IniReadString ('stop', festIni, section = 'time')
ELSE
    CALL Catch ('error', 'Fest', 'stop missing in time' )
END IF

timeStop = timeString

CALL Catch ('info', 'Fest', 'stop time: ', argument = timeString )

!------------------------------------------------------------------------------
!                           result folder
!------------------------------------------------------------------------------
CALL Catch ('info', 'Fest', 'result folder configuration' )

IF ( verbose ) THEN
    WRITE(*,*)
    WRITE(*,*) '...setting folder for results...'
END IF

!check whether section is present
IF ( .NOT. SectionIsPresent ('result', festIni) ) THEN  !section is mandatory
    CALL Catch ('error', 'Fest', 'result section missing in configuration file' )
END IF    

!folder
IF (KeyIsPresent('folder', festIni, section = 'result')) THEN	
   path_out = IniReadString ('folder', festIni, section = 'result')
ELSE
    CALL Catch ('error', 'Fest', 'folder missing in result' )
END IF

CALL Catch ('info', 'Fest', 'result folder: ', argument = TRIM(path_out) )


!------------------------------------------------------------------------------
!                        configure simulation domain 
!------------------------------------------------------------------------------
CALL Catch ('info', 'Fest', 'simulation domain configuration' )

IF ( verbose ) THEN
    WRITE(*,*)
    WRITE(*,*) '...configuring simulation domain...'
END IF
    
!check whether section is present
IF ( .NOT. SectionIsPresent ('domain', festIni) ) THEN !section is mandatory
    CALL Catch ('error', 'Fest', 'domain section missing in configuration file' )
END IF  

!find configuration file
IF (KeyIsPresent('conf-file', festIni, section = 'domain')) THEN	
   string = IniReadString ('conf-file', festIni, section = 'domain')
ELSE
    CALL Catch ('error', 'Fest', 'conf-file missing in domain' )
END IF

!read data
CALL DomainInit (string)

!------------------------------------------------------------------------------
!                        morphological properties
!------------------------------------------------------------------------------ 

 IF ( verbose ) THEN
        WRITE(*,*)
        WRITE(*,*) '...configuring morphological properties...'
    END IF

!check whether section is present. Optional section
IF ( .NOT. SectionIsPresent ('morphology', festIni) ) THEN !section is mandatory
    CALL Catch ('info', 'Fest', 'morphological properties configuration' ) 
END IF  

!find configuration file
IF (KeyIsPresent('conf-file', festIni, section = 'morphology')) THEN	
   string = IniReadString ('conf-file', festIni, section = 'morphology')
   !read data
   CALL MorphologyInit (string, mask)
ELSE
    CALL Catch ('error', 'Fest', 'conf-file missing in morphology' )
END IF


!------------------------------------------------------------------------------
!                  configure meteo input
!------------------------------------------------------------------------------
CALL Catch ('info', 'Fest', 'meteo input configuration' ) 

IF ( verbose ) THEN
    WRITE (*,*)
    WRITE (*,*) '...configuring meteorological forcings...'
END IF

!check whether section is present
IF ( .NOT. SectionIsPresent ('meteo', festIni) ) THEN !section is mandatory
    CALL Catch ('error', 'Fest', 'meteo section missing in configuration file')
END IF  

!time step
IF (KeyIsPresent('dt', festIni, section = 'meteo')) THEN	
   dtMeteo = IniReadInt ('dt', festIni, section = 'meteo')
ELSE
    CALL Catch ('error', 'Fest', 'dt missing in meteo' )
END IF

!find configuration file
IF (KeyIsPresent('conf-file', festIni, section = 'meteo')) THEN	
   string = IniReadString ('conf-file', festIni, section = 'meteo')
ELSE
    CALL Catch ('error', 'Fest', 'conf-file missing in meteo' )
END IF

!read data
CALL MeteoInit (string, timeStart, mask, dem, dem_loaded, albedo_loaded)


!spatial data out dt
IF (KeyIsPresent('dt-out-spatial', festIni, section = 'meteo')) THEN	
   dtOutMeteo = IniReadInt ('dt-out-spatial', festIni, section = 'meteo')
   IF ( dtOutMeteo > 0 ) THEN
       timeOutMeteo = timeStart
   ELSE
       timeOutMeteo = timeStart + (-1)
   END IF
ELSE
    CALL Catch ('warning', 'Fest', 'dt-out-spatial missing in meteo' )
    dtOutMeteo = - 1
    timeOutMeteo = timeStart + (-1)
END IF

!point data out
IF (KeyIsPresent('out-point-file', festIni, section = 'meteo')) THEN	
   string = IniReadString ('out-point-file', festIni, section = 'meteo')
   CALL MeteoPointInit (string, path_out, timeStart)
END IF	

timeNewMeteo = timeStart


!------------------------------------------------------------------------------
!                         irrigation 
!------------------------------------------------------------------------------
CALL Catch ('info', 'Fest', 'irrigation configuration' ) 

!check whether section is present
IF ( SectionIsPresent ('irrigation', festIni) ) THEN !section is optional
    
    IF ( verbose ) THEN
      WRITE(*,*)
      WRITE(*,*) '...configuring irrigation...'
    END IF
    
    !time step
    IF (KeyIsPresent('dt', festIni, section = 'irrigation')) THEN	
       dtIrrigation = IniReadInt ('dt', festIni, section = 'irrigation')
    ELSE
        CALL Catch ('error', 'Fest', 'dt missing in irrigation' )
    END IF
    
    !dt output
    IF (KeyIsPresent('dt-out', festIni, section = 'irrigation')) THEN	
       dtOutIrrigation = IniReadInt ('dt-out', festIni, section = 'irrigation')
       IF ( dtOutIrrigation > 0 ) THEN
          timeOutIrrigation = timeStart
       ELSE
          timeOutIrrigation = timeStart + (-1)
       END IF
    ELSE
        CALL Catch ('error', 'Fest', 'dt-out missing in irrigation' )
    END IF
    
    !configuration file
    IF ( dtIrrigation > 0 ) THEN
        IF (KeyIsPresent('conf-file', festIni, section = 'irrigation')) THEN	
           string = IniReadString ('conf-file', festIni, section = 'irrigation')
           CALL IrrigationConfig (string, path_out, dtOutIrrigation )
        ELSE
            CALL Catch ('error', 'Fest', 'conf-file missing in irrigation' )
        END IF
        
        timeNewIrrigation = timeStart
        
    END IF
ELSE
    dtIrrigation = 0
END IF  
	
!------------------------------------------------------------------------------
!                         canopy interception 
!------------------------------------------------------------------------------

!_________________________________________
!  da implementare che se c'è neve non calcolo intercettazione BO! serve?
!________________________


CALL Catch ('info', 'Fest', 'canopy interception configuration' ) 

!raineff is used even when canopy interception is not simulated
CALL NewGrid (raineff, mask, 0.)

!check whether section is present
IF ( SectionIsPresent ('canopy-interception', festIni) ) THEN !section is optional
    
    IF ( verbose ) THEN
      WRITE(*,*)
      WRITE(*,*) '...configuring canopy interception...'
    END IF
    
    !dt
    IF (KeyIsPresent('dt', festIni, section = 'canopy-interception')) THEN	
       dtCanopyInterception = IniReadInt ('dt', festIni, section = 'canopy-interception')
    ELSE
        CALL Catch ('error', 'Fest', 'dt missing in canopy-interception' )
    END IF
    
    IF ( dtCanopyInterception > 0 ) THEN
        !configuration file
        IF (KeyIsPresent('conf-file', festIni, section = 'canopy-interception')) THEN	
           string = IniReadString ('conf-file', festIni, section = 'canopy-interception')
           CALL InterceptionInit ( string, mask )
        ELSE
            CALL Catch ('error', 'Fest', 'conf-file missing in canopy-interception' )
        END IF
    
        !spatial data out dt
        IF (KeyIsPresent('dt-out-spatial', festIni, section = 'canopy-interception')) THEN	
           dtOutCanopy = IniReadInt ('dt-out-spatial', festIni, section = 'canopy-interception')
        ELSE
            CALL Catch ('error', 'Fest', 'dt-out-spatial missing in canopy-interception' )
        END IF
        
        timeNewCanopyInterception = timeStart
    
        IF (dtOutCanopy > 0) THEN
            timeOutCanopy = timeStart
        ELSE
             timeOutCanopy = timeStart + (-1)
        END IF
    ELSE
        dtOutCanopy = 0
        timeNewCanopyInterception = timeStart + (-1)
        timeOutCanopy = timeStart + (-1)
    END IF
    
ELSE
   dtCanopyInterception = 0 
   dtOutCanopy = 0
   timeNewCanopyInterception = timeStart + (-1)
   timeOutCanopy = timeStart + (-1)
END IF  
	

!------------------------------------------------------------------------------
!                 Plants dynamic simulation
!------------------------------------------------------------------------------
CALL Catch ('info', 'Fest', 'plants configuration' ) 

!check whether section is present
IF ( SectionIsPresent ('plants', festIni) ) THEN !section is optional
    
    IF ( verbose ) THEN
      WRITE(*,*)
      WRITE(*,*) '...configuring plants simulation...'
   END IF
    
    !dt
    IF (KeyIsPresent('dt', festIni, section = 'plants')) THEN	
       dtPlants = IniReadInt ('dt', festIni, section = 'plants')
    ELSE
        CALL Catch ('error', 'Fest', 'dt missing in plants' )
    END IF
    
     !configuration file
    IF (KeyIsPresent('conf-file', festIni, section = 'plants')) THEN	
        string = IniReadString ('conf-file', festIni, section = 'plants')
        CALL PlantsConfig (string, mask, timeStart, timeStop)
    ELSE
        CALL Catch ('error', 'Fest', 'conf-file missing in plants' )
    END IF
    
    !spatial data out dt
    IF (KeyIsPresent('dt-out-spatial', festIni, section = 'plants')) THEN	
        dtOutPlants = IniReadInt ('dt-out-spatial', festIni, section = 'plants')
    ELSE
        CALL Catch ('warning', 'Fest', 'dt-out-spatial missing in plants' )
        dtOutPlants = 0
    END IF
    
    IF (dtOutPlants > 0) THEN
        timeOutPlants = timeStart
    END IF
    
    IF ( dtPlants > 0 ) THEN
        timeNewPlants = timeStart
    ELSE
        dtOutPlants = 0
        timeNewPlants = timeStart + (-1)
        timeOutPlants = timeStart + (-1)
    END IF

ELSE
    dtPlants = 0
    dtOutPlants = 0
    timeNewPlants = timeStart + (-1)
    timeOutPlants = timeStart + (-1)
END IF  
	

!------------------------------------------------------------------------------
!               snow accumulation and melting
!------------------------------------------------------------------------------
CALL Catch ('info', 'Fest', 'snow configuration' ) 

!rainfallRate is used even when snow is not simulated
CALL NewGrid (rainfallRate, mask, 0.)

!check whether section is present
IF ( SectionIsPresent ('snow', festIni) ) THEN !section is optional
    
    IF ( verbose ) THEN
        WRITE(*,*)
        WRITE(*,*) '...configuring snow accumulation and melting...'
    END IF
    
    !dt
    IF (KeyIsPresent('dt', festIni, section = 'snow')) THEN	
       dtSnow = IniReadInt ('dt', festIni, section = 'snow')
    ELSE
        CALL Catch ('error', 'Fest', 'dt missing in snow' )
    END IF
    
    
    IF (dtSnow > 0) THEN
        
		!check dt
		IF ( dtSnow /= dtMeteo ) THEN
			dtSnow=dtMeteo
			WRITE(*,*)'WARNING: dtSnow \E8 posto uguale a dt_acquisizione_meteo'
		END IF
		timeNewSnow = timeStart
        
        !configuration file
        IF (KeyIsPresent('conf-file', festIni, section = 'snow')) THEN	
        string = IniReadString ('conf-file', festIni, section = 'snow')
        CALL SnowInit (string, mask, timeStart )
        ELSE
        CALL Catch ('error', 'Fest', 'conf-file missing in snow' )
        END IF
        
        !spatial data out dt
        IF (KeyIsPresent('dt-out-spatial', festIni, section = 'snow')) THEN	
           dtOutSnow = IniReadInt ('dt-out-spatial', festIni, section = 'snow')
        ELSE
            CALL Catch ('error', 'Fest', 'dt-out-spatial missing in snow' )
        END IF
        
        IF ( dtOutSnow > 0 ) THEN
			timeOutSnow = timeStart
        END IF
        
        !point data out
		IF (KeyIsPresent('out-point-file', festIni, section = 'snow')) THEN	
		   string = IniReadString ('out-point-file', festIni, section = 'snow')
		   CALL SnowPointInit ( string, path_out, timeStart )
		END IF
        
    ELSE
        dtOutSnow = 0
        timenEWSnow = timeStart + (-1)
        timeOutSnow = timeStart + (-1)
    END IF
ELSE    
    dtSnow = 0
    dtOutSnow = 0
    timenEWSnow = timeStart + (-1)
    timeOutSnow = timeStart + (-1)
END IF  


!------------------------------------------------------------------------------
!                      glacier
!------------------------------------------------------------------------------
CALL Catch ('info', 'Fest', 'glacer configuration' ) 

!check whether section is present
IF ( SectionIsPresent ('glacier', festIni) ) THEN !section is optional
    
    IF ( verbose ) THEN
        WRITE(*,*)
        WRITE(*,*) '...configuring glacier accumulation and ablation...'
    END IF
    
    !dt
    IF (KeyIsPresent('dt', festIni, section = 'glacier')) THEN	
       dtIce = IniReadInt ('dt', festIni, section = 'glacier')
    ELSE
        CALL Catch ('error', 'Fest', 'dt missing in glacier' )
    END IF
    
    
    IF (dtIce > 0) THEN
        
		!check dt
		IF ( dtIce /= dtMeteo ) THEN
			dtIce = dtMeteo
            CALL Catch ('info', 'Fest', 'dtIce set equal to dtMeteo' )
        END IF
        
		timeNewIce = timeStart
        
        !configuration file
        IF (KeyIsPresent('conf-file', festIni, section = 'glacier')) THEN	
           string = IniReadString ('conf-file', festIni, section = 'glacier')
           CALL GlacierInit (string, mask, timeStart )
        ELSE
           CALL Catch ('error', 'Fest', 'conf-file missing in glacier' )
        END IF
        
        !spatial data out dt
        IF (KeyIsPresent('dt-out-spatial', festIni, section = 'glacier')) THEN	
           dtOutIce = IniReadInt ('dt-out-spatial', festIni, section = 'glacier')
        ELSE
            CALL Catch ('error', 'Fest', 'dt-out-spatial missing in glacier' )
        END IF
        
        IF ( dtOutIce > 0 ) THEN
			timeOutIce = timeStart
        END IF
        
        !point data out
		IF (KeyIsPresent('out-point-file', festIni, section = 'glacier')) THEN	
		   string = IniReadString ('out-point-file', festIni, section = 'glacier')
           CALL GlacierPointInit ( string, path_out, timeStart )
		END IF
        
    ELSE
         dtOutIce = 0
         timeNewIce = timeStart + (-1)
         timeOutIce = timeStart + (-1)
	END IF
ELSE
    dtIce = 0
    dtOutIce = 0
    timeNewIce = timeStart + (-1)
    timeOutIce = timeStart + (-1)
END IF  
 
!------------------------------------------------------------------------------
!                              Soil balance
!------------------------------------------------------------------------------
CALL Catch ('info', 'Fest', 'soil balance configuration' ) 


!check whether section is present
IF ( SectionIsPresent ('soil-balance', festIni) ) THEN !section is optional
    
    IF ( verbose ) THEN
       WRITE(*,*)
       WRITE(*,*) '...configuring soil water balance...'
    END IF
    
    !dt
    IF (KeyIsPresent('dt', festIni, section = 'soil-balance')) THEN	
       dtSoilBalance = IniReadInt ('dt', festIni, section = 'soil-balance')
    ELSE
        CALL Catch ('error', 'Fest', 'dt missing in soil-balance' )
    END IF
    
    IF ( dtSoilBalance  > 0 ) THEN
        
        timeNewSoilBalance = timeStart
        
        !configuration file
        IF (KeyIsPresent('conf-file', festIni, section = 'soil-balance')) THEN	
           string = IniReadString ('conf-file', festIni, section = 'soil-balance')
           CALL InitSoilBalance (string, flowDirection, timeStart)
        ELSE
            CALL Catch ('error', 'Fest', 'conf-file missing in soil-balance' )
        END IF
    
        !spatial data out dt
        IF (KeyIsPresent('dt-out-spatial', festIni, section = 'soil-balance')) THEN	
           dtOutSoilBalance = IniReadInt ('dt-out-spatial', festIni, &
                                           section = 'soil-balance')
        ELSE
            dtOutSoilBalance = 0
        END IF
    
        IF (dtOutSoilBalance > 0) THEN
            timeOutSoilBalance = timeStart
        END IF
    ELSE
        dtOutSoilBalance = 0
        timeNewSoilBalance = timeStart + (-1)
        timeOutSoilBalance = timeStart + (-1)
    END IF
ELSE
    dtSoilBalance = 0
    dtOutSoilBalance = 0
    timeNewSoilBalance = timeStart + (-1)
    timeOutSoilBalance = timeStart + (-1)
END IF  
    
!------------------------------------------------------------------------------
!                            groundwater
!------------------------------------------------------------------------------
CALL Catch ('info', 'Fest', 'groundwater configuration' )

!check whether section is present
IF ( SectionIsPresent ('groundwater', festIni) ) THEN !section is optional
    
    IF ( verbose ) THEN
      WRITE(*,*)
      WRITE(*,*) '...configuring groundwater...'
    END IF
    
    !dt
    IF (KeyIsPresent('dt', festIni, section = 'groundwater')) THEN	
       dtGroundwater = IniReadInt ('dt', festIni, section = 'groundwater')
    ELSE
        CALL Catch ('error', 'Fest', 'dt missing in groundwater' )
    END IF
    
    IF ( dtGroundwater > 0 ) THEN
        
         timeNewGroundwater = timeStart	
        
        !configuration file
        IF (KeyIsPresent('conf-file', festIni, section = 'groundwater')) THEN	
           string = IniReadString ('conf-file', festIni, section = 'groundwater')
       
           CALL GroundwaterInit (string)
      				
        ELSE
            CALL Catch ('error', 'Fest', 'conf-file missing in groundwater' )
        END IF
    
        !spatial data out dt
        IF (KeyIsPresent('dt-out-aquifer', festIni, section = 'groundwater')) THEN	
           dtOutGroundwater = IniReadInt ('dt-out-aquifer', festIni, section = 'groundwater')
           IF ( dtOutGroundwater > 0 ) THEN
                timeOutGroundwater = timeStart
                CALL GroundwaterOutputInit ( path_out)
           ELSE
               timeOutGroundwater = timeStart + (-1)
           END IF
        ELSE
            CALL Catch ('error', 'Fest', 'dt-out-aquifer missing in groundwater' )
        END IF
    
        !point data out 
        IF (KeyIsPresent('out-point-file', festIni, section = 'groundwater')) THEN	
            string = IniReadString ('out-point-file', festIni, section = 'groundwater')
    
            CALL GroundwaterPointInit ( string, path_out, timestart )

        END IF
    
    ELSE
         dtOutGroundwater = 0
         timeNewGroundwater = timeStart + (-1)
         timeOutGroundwater = timeStart + (-1)
    END IF
ELSE
    dtGroundwater = 0
    dtOutGroundwater = 0
    timeNewGroundwater = timeStart + (-1)
    timeOutGroundwater = timeStart + (-1)
END IF  	


!------------------------------------------------------------------------------
!                   soil erosion and sediment transport
!------------------------------------------------------------------------------

CALL Catch ('info', 'Fest', 'sediment erosion configuration' ) 

!check whether section is present
IF ( SectionIsPresent ('sediment', festIni) ) THEN !section is optional
    
    IF ( verbose ) THEN
      WRITE(*,*)
      WRITE(*,*) '...configuring soil erosion and sediment transport...'
    END IF
    
    !dt
    IF (KeyIsPresent('dt', festIni, section = 'sediment')) THEN	
       dtCalcSediment = IniReadInt ('dt', festIni, section = 'sediment')
    ELSE
        CALL Catch ('error', 'Fest', 'dt missing in sediment' )
    END IF
    
    IF ( dtCalcSediment > 0 ) THEN
    
        !configuration file
        IF (KeyIsPresent('conf-file', festIni, section = 'sediment')) THEN	
           string = IniReadString ('conf-file', festIni, section = 'sediment')
           CALL SedimentInit (string, dtCalcSedimentRouting, string)
           timeUpdateSediment = timeStart
           timeUpdateSedimentRouting = timeStart
        ELSE
            CALL Catch ('error', 'Fest', 'conf-file missing in sediment' )
        END IF
    
        !spatial data out dt
        IF (KeyIsPresent('dt-out-spatial', festIni, section = 'sediment')) THEN	
           dtOutSediment = IniReadInt ('dt-out-spatial', festIni, section = 'sediment')
        ELSE
            CALL Catch ('error', 'Fest', 'dt-out-spatial missing in sediment' )
        END IF
    
        IF (dtOutSediment > 0) THEN
            timeOutSoilBalance = timeStart
        END IF
    
        timeOutSediment = timeStart
    ELSE
        dtOutSediment = 0
        timeUpdateSediment = timeStart + (-1)
        timeOutSediment = timeStart + (-1)
    END IF
ELSE
    dtCalcSediment = 0
    dtOutSediment = 0
    timeUpdateSediment = timeStart + (-1)
    timeOutSediment = timeStart + (-1)
END IF
    
    !TODO
    
  !      IF (routeSediment) THEN
		!	!set file for output of sediment routing
		!	fileunit_out_sediment_routing = GetUnit ()
		!	OPEN (fileunit_out_sediment_routing,file=string)
		!	WRITE(*,*)'nome file anagrafica out punti propagazione sedimenti: ',string(1:LEN_TRIM(string))
		!	CALL leggi_anagrafica(xsOutSedimentRouting,fileunit_out_sediment_routing)
		!	!check dt out sediment routing
		!	IF(.not.(mod(xsOutSedimentRouting%cadenza_misure,dtCalcSedimentRouting)==0)) THEN
		!		xsOutSedimentRouting%cadenza_misure=dtCalcSedimentRouting
		!		WRITE(*,*)'WARNING: xsOutSedimentRouting%cadenza_misure non multiplo di dtCalcSedimentRouting'
		!		WRITE(*,*)'         si impone: xsOutSedimentRouting%cadenza_misure = dtCalcSedimentRouting'
		!	END IF
		!	WRITE(*,*)xsOutSedimentRouting%num_staz,' sezioni per risultati propagazione sedimenti'
		!	DO k=1,xsOutSedimentRouting%num_staz
		!		xsOutSedimentRouting%staz_oss(k)%misura=0.
		!	END DO
		!	CALL righe_colonne(xsOutSedimentRouting,grid_r=QoutSed)
		!	CLOSE (fileunit_out_sediment_routing)
  !  			
		!	!write header in output file
		!	fileunit_out_sediment_routing = GetUnit ()
		!	OPEN (fileunit_out_sediment_routing,file=path_out(1:LEN_TRIM(path_out))//'sediment_routing.fest')
		!	WRITE(fileunit_out_sediment_routing,'(a)')'FEST-EWB: stima dei sedimentogrammi'
		!	WRITE(fileunit_out_sediment_routing,'(a)')'anagrafica'
		!	CALL scrivi_anagrafica(xsOutSedimentRouting,fileunit_out_sediment_routing)
		!	WRITE(fileunit_out_sediment_routing,'(a)')
		!	WRITE(fileunit_out_sediment_routing,'(a)')
		!	WRITE(fileunit_out_sediment_routing,'(a)')'dati'
		!	WRITE(fileunit_out_sediment_routing,'(a,"	",\)')'#'
		!	DO k=1,xsOutSedimentRouting%num_staz-1
		!		WRITE(fileunit_out_sediment_routing,'(a,"	",\)')(xsOutSedimentRouting%staz_oss(k)%denominazione)
		!	END DO
		!	WRITE(fileunit_out_sediment_routing,'(a)')(xsOutSedimentRouting%staz_oss(xsOutSedimentRouting%num_staz)%denominazione)
  !  			
		!	!time
		!	timeOutSedimentRouting = timeStart
		!END IF
  !  
    
    

    
!------------------------------------------------------------------------------
!                            spatial average
!------------------------------------------------------------------------------

CALL Catch ('info', 'Fest', 'spatial-average configuration' ) 


!check whether section is present
IF ( SectionIsPresent ('spatial-average', festIni) ) THEN !section is optional
    
IF ( verbose ) THEN
    WRITE(*,*)
    WRITE(*,*) '...configuring spatial average for output...'
END IF
    
!configuration file
IF (KeyIsPresent('conf-file', festIni, section = 'spatial-average')) THEN	
    string = IniReadString ('conf-file', festIni, section = 'spatial-average')
       

    CALL InitSpatialAverageMeteo (string, path_out, temperatureAir, &
                        temperatureAirDailyMean, temperatureAirDailyMax, &
                        temperatureAirDailyMin, precipitationRate, &
                        relHumidityAir, radiation,netRadiation, &
                        windSpeed, precipitationRateDaily, irrigationFlux ) 


    CALL InitSpatialAverageBalance (string, path_out, soilMoisture, &
                                    soilMoistureRZ, soilMoistureRZ, &
                                    runoff, infilt, percolation, &
                                    et, pet, caprise, balanceError)

    CALL InitSpatialAverageSnow   (string, path_out, rainfallRate, swe, &
                                    meltCoefficient, waterInSnow, snowMelt ) 

    CALL InitSpatialAverageGlaciers   (string, path_out, icewe, waterInIce, iceMelt )  

    CALL InitSpatialAverageSediment   (string, path_out, interrillErosion) 
        
    CALL InitSpatialAverageCanopy   (string, path_out, canopyStorage, raineff, canopyPT ) 
        
    CALL InitSpatialAveragePlants   (string, path_out, lai, gpp, npp, carbonstem, &
                                    carbonroot, carbonleaf, fvcover, dbh, plantsHeight, &
                                    density, stemyield )
ELSE
    CALL Catch ('error', 'Fest', 'conf-file missing in spatial-average' )
END IF
    
!initialize  meteo output file
IF ( dtOutMeteo > 0 ) THEN
                   
    CALL ComputeSpatialAverageMeteo (dtMeteo, temperatureAir, &
                            temperatureAirDailyMean, &
                            temperatureAirDailyMax, temperatureAirDailyMin, &
                            precipitationRate, relHumidityAir, radiation, &
                            netRadiation, windSpeed, precipitationRateDaily, &
                            irrigationFlux)
				
    CALL ExportSpatialAverageMeteo ( init = .TRUE. )
END IF
    
!initialize snow output file
IF ( dtOutSnow > 0 ) THEN
        
    CALL ComputeSpatialAverageSnow  ( dtOutSnow, rainfallRate, &
                                    swe, meltCoefficient, waterInSnow, &
                                    snowMelt )
        
    CALL ExportSpatialAverageSnow (init = .TRUE.)		
END IF
    
!initialize glaciers output file
IF ( dtOutIce > 0 ) THEN
    CALL ComputeSpatialAverageGlaciers  (dtOutIce, icewe, &
                                            waterInIce, iceMelt ) 
    CALL ExportSpatialAverageGlaciers (init = .TRUE.)
END IF
    
!initailize vegetation canopy output file
IF (dtOutCanopy > 0) THEN
    CALL ExportSpatialAverageCanopy (init = .TRUE.)
END IF
    
!initailize plants output file
IF (dtOutPlants > 0) THEN
    CALL ExportSpatialAveragePlants (init = .TRUE.)
END IF
    
!initailize soil balance output file
IF (dtOutSoilBalance>0) THEN	
    CALL ExportSpatialAverageBalance (init = .TRUE.)
END IF
    
!initailize sediment output file
IF (dtCalcSediment > 0) THEN	
    CALL ExportSpatialAverageSediment (init = .TRUE.)
END IF
    
END IF  




!------------------------------------------------------------------------------
!                            basin average
!------------------------------------------------------------------------------

CALL Catch ('info', 'Fest', 'basin-average configuration' ) 


!check whether section is present
IF ( SectionIsPresent ('basin-average', festIni) ) THEN !section is optional
    
    IF ( verbose ) THEN
        WRITE(*,*)
        WRITE(*,*) '...configuring basin average for output...'
    END IF

    !output  file
    IF (KeyIsPresent('out-point-file', festIni, section = 'basin-average')) THEN	
        string = IniReadString ('out-point-file', festIni, section = 'basin-average')
        CALL ReadPointFileBasinAverage ( string )
    ELSE
        CALL Catch ('error', 'Fest', 'out-point-file missing in basin-average' )
    END IF
    
    !configuration file
    IF (KeyIsPresent('conf-file', festIni, section = 'basin-average')) THEN	
        string = IniReadString ('conf-file', festIni, section = 'basin-average')
       
        CALL InitBasinAverage (string, path_out, temperatureAir, &
                            temperatureAirDailyMean, temperatureAirDailyMax, &
                            temperatureAirDailyMin, precipitationRate, &
                            relHumidityAir, radiation,netRadiation, &
                            windSpeed, precipitationRateDaily, irrigationFlux, &
                            swe, soilMoisture, runoff, infilt, percolation, &
                            et, pet, deltaSoilMoisture, snowMelt, icewe, &
                            iceMelt) 
    
    !set first export time
    timeNewBasinAverage = timeStart
    
    ELSE
        CALL Catch ('error', 'Fest', 'conf-file missing in basin-average' )
    END IF

ELSE

    timeNewBasinAverage = timeStart + (-1)
    
END IF


!------------------------------------------------------------------------------
!                            raster export
!------------------------------------------------------------------------------

CALL Catch ('info', 'FeST', 'raster-export configuration' ) 


!check whether section is present
IF ( SectionIsPresent ('raster-export', festIni) ) THEN !section is optional
    
    IF ( verbose ) THEN
        WRITE(*,*)
        WRITE(*,*) '...configuring raster export for output...'
    END IF
    
    !configuration file
    IF (KeyIsPresent('conf-file', festIni, section = 'raster-export')) THEN	
        string = IniReadString ('conf-file', festIni, section = 'raster-export')
       
        CALL InitRasterExport (string, temperatureAir, precipitationRate, &
                            relHumidityAir, radiation, netRadiation, &
                            windSpeed,swe, soilMoisture, runoff, infilt, &
                            percolation, et, pet ) 
    
    !set first time
    timeNewRasterExport = timeStart
    
    ELSE
        CALL Catch ('error', 'Fest', 'conf-file missing in raster-export' )
    END IF

ELSE

    timeNewRasterExport = timeStart + (-1)
    
END IF






!------------------------------------------------------------------------------
!                            discharge routing
!------------------------------------------------------------------------------

CALL Catch ('info', 'Fest', 'surface routing configuration' ) 


!check whether section is present
IF ( SectionIsPresent ('discharge-routing', festIni) ) THEN !section is optional
    
    IF ( verbose ) THEN
      WRITE(*,*)
      WRITE(*,*) '...configuring discharge routing...'
   END IF
    
    !dt
    IF (KeyIsPresent('dt', festIni, section = 'discharge-routing')) THEN	
       dtDischargeRouting = IniReadInt ('dt', festIni, section = 'discharge-routing')
    ELSE
        CALL Catch ('error', 'Fest', 'dt missing in discharge-routing' )
    END IF
    
    IF ( dtDischargeRouting > 0 ) THEN
        
        timeNewDischargeRouting = timeStart
     
        !configuration file
        IF (KeyIsPresent('conf-file', festIni, section = 'discharge-routing')) THEN	
           string = IniReadString ('conf-file', festIni, section = 'discharge-routing')
           CALL  InitDischargeRouting (fileini = string, time = timeStart, &
                            path_out = path_out, &
                        flowdirection = flowDirection, &
                        domain = mask , &
                        storage = storageChannelSurf, &
                        netRainfall = runoff, &
                        dtrk = dtReservoir )		
        ELSE
            CALL Catch ('error', 'Fest', 'conf-file missing in discharge-routing' )
        END IF
    
        !output data file
        IF (KeyIsPresent('out-point-file', festIni, section = 'discharge-routing')) THEN	
            string = IniReadString ('out-point-file', festIni, section = 'discharge-routing')
            CALL DischargePointInit (string, path_out, timeStart)
        END IF
    ELSE
        timeNewDischargeRouting = timeStart + (-1)
    END IF
    
ELSE
    dtDischargeRouting = 0
    timeNewDischargeRouting = timeStart + (-1)
END IF 

!------------------------------------------------------------------------------
!                             hot start
!------------------------------------------------------------------------------
CALL Catch ('info', 'Fest', 'save-hot-start configuration' ) 

!check whether section is present
IF ( SectionIsPresent ('save-hot-start', festIni) ) THEN !section is optional
    
     hotstart = .TRUE.
    
    IF ( verbose ) THEN
      WRITE(*,*)
      WRITE(*,*) '...configuring hot start...'
    END IF
    
    !time 
    IF (KeyIsPresent('time', festIni, section = 'save-hot-start')) THEN	
       cron = IniReadString ('time', festIni, section = 'save-hot-start')
       
       CALL CronParseString (cron, cronHotstart )
       
       !CALL ConfigureHotStart (cron)
       
    ELSE
        CALL Catch ('error', 'Fest', 'time missing in save-hot-start' )
    END IF
    
    !folder
    IF (KeyIsPresent('folder', festIni, section = 'save-hot-start')) THEN	
        pathHotstart = IniReadString ('folder', festIni, section = 'save-hot-start')
    ELSE
        CALL Catch ('error', 'Fest', 'folder missing in save-hot-start' )
    END IF
    
    
    !save last time 
    IF (KeyIsPresent('save-last', festIni, section = 'save-hot-start')) THEN	
       IF ( IniReadInt ('save-last', festIni, section = 'save-hot-start') == 1 ) THEN
           saveLast = .TRUE.
       ELSE
           saveLast = .FALSE.
       END IF
    ELSE !set to default
        saveLast = .FALSE.
    END IF
    
    !IF ( dtHotstart > 0 ) THEN
    !    
    !    timeNewHotstart = timeStart
    !    
    !    !folder
    !    IF (KeyIsPresent('folder', festIni, section = 'save-hot-start')) THEN	
    !       pathHotstart = IniReadString ('folder', festIni, section = 'save-hot-start')
    !    ELSE
    !        CALL Catch ('error', 'Fest', 'folder missing in save-hot-start' )
    !    END IF
    !
    !ELSE
    !    timeNewHotstart = timeStart + (-1)
    !END IF
!ELSE
!    dtHotstart = 0
!    timeNewHotstart = timeStart + (-1)
    
ELSE
     hotstart = .FALSE.
END IF

IF ( verbose ) THEN
  WRITE(*,*)
  WRITE(*,*) '....configuration completed!'
END IF


!------------------------------------------------------------------------------
!                    set calculation time step
!------------------------------------------------------------------------------

dtArray (1) = dtMeteo 
dtArray (2) = dtSoilBalance
dtArray (3) = dtGroundwater
dtArray (4) = dtDischargeRouting

dtCalc = HUGE ( dtCalc )
DO i = 1, SIZE (dtArray)
   IF ( dtArray (i) > 0 .AND. dtArray (i) < dtCalc ) THEN
       dtCalc = dtArray (i)
   END IF
END DO


!------------------------------------------------------------------------------
!                 start simulation
!------------------------------------------------------------------------------

IF ( verbose ) THEN
    WRITE(*,*)
    WRITE(*,*) 'starting simulation...'
    WRITE(*,*)
    WRITE(*,'(a,I10,a)') 'computation time step: ', dtCalc, '  second'
END IF

time = timeStart

!DO WHILE( .NOT. ( time < timeStop) )
DO WHILE(  time <= timeStop )
	timeString = time
    
	IF ( verbose ) THEN    
		WRITE(*,'(a25)') timeString
	END IF

	
	!-----------------------------------------------------------------------------
	!   read meteorological forcings
	!-----------------------------------------------------------------------------
	IF (time == timeNewMeteo) THEN
        IF ( verbose ) THEN 
			WRITE(*,*)'read meteorological forcings'
        END IF
      
        CALL MeteoRead (time,dem, albedo)
      
            
		timeNewMeteo = timeNewMeteo + dtMeteo
                 
        IF ( dtOutMeteo > 0 ) THEN
            timeSpatialAverageMeteo = time 
            CALL ComputeSpatialAverageMeteo (dtMeteo, temperatureAir, &
                                        temperatureAirDailyMean, &
                                        temperatureAirDailyMax, temperatureAirDailyMin, &
                                        precipitationRate, relHumidityAir, radiation, &
                                        netRadiation, windSpeed, precipitationRateDaily, &
                                        irrigationFlux)
        END IF
			
    END IF
    
    !-----------------------------------------------------------------------------
	!   update irrigation
	!-----------------------------------------------------------------------------
    IF (time == timeNewIrrigation) THEN

        IF ( verbose ) THEN 
            WRITE(*,*) 'update irrigation'
        END IF
        
		CALL IrrigationUpdate (time, soilSat, Qin, dtOutIrrigation, &
                               timeOutIrrigation)
        
        timeNewIrrigation = timeNewIrrigation + dtIrrigation
    END IF  
        
	!-----------------------------------------------------------------------------
	!                          snow
	!-----------------------------------------------------------------------------
	IF (dtSnow>0) THEN
		IF (time == timeNewSnow) THEN
            IF ( verbose ) THEN 
			   WRITE(*,*)'snow'
            END IF
            
            CALL SnowUpdate (time, mask)
           
			timeNewSnow = timeNewSnow + dtSnow
           
            timeSpatialAverageSnow = time
            CALL ComputeSpatialAverageSnow  (dtOutSnow, rainfallRate, swe, &
                                       meltCoefficient, waterInSnow, &
                                       snowMelt)
		END IF
    ELSE
        
        rainfallRate = precipitationRate
		
    END  IF

	!-----------------------------------------------------------------------------
	!                      glacier
	!-----------------------------------------------------------------------------
    
	IF ( time == timeNewIce ) THEN
        IF ( verbose ) THEN 
		    WRITE(*,*) 'glaciers'
        END IF
        
        CALL GlacierUpdate (time, mask, rainfallRate ) 
		
		timeNewIce = timeNewIce + dtIce
    
        timeSpatialAverageIce = time
        CALL ComputeSpatialAverageGlaciers  (dtOutIce, icewe, &
                                             waterInIce, iceMelt) 
    END IF
    
	!-----------------------------------------------------------------------------
	!CANOPY INTERCEPTION
	!-----------------------------------------------------------------------------
    
    IF ( dtCanopyInterception > 0 ) THEN
        
        IF ( time == timeNewCanopyInterception) THEN
            IF ( verbose ) THEN 
                WRITE(*,*) 'canopy interception'
            END IF
        
            CALL Throughfall (rainfallRate , lai, mask, fvcover, raineff)
        
            !CALL AdjustPT(fvcover, mask, tp_oss, dtetp)
        
            timeNewCanopyInterception = timeNewCanopyInterception + dtCanopyInterception
            
            timeSpatialAverageCanopy = time
            CALL ComputeSpatialAverageCanopy  (dtCanopyInterception, canopyStorage, raineff, canopyPT)
            
        END IF
        
    ELSE !copy rainfallRate into raineff
        
        raineff = rainfallRate
        
    END IF
    
    
	!-------------------------------------------------------------------------------------------------------------
	!FOREST DYNAMIC
	!-------------------------------------------------------------------------------------------------------------
    !check if plants parameters need update
    IF ( updatePlantsParameters ) THEN
        CALL PlantsParameterUpdate ( time )
    END IF
    
    
     IF ( dtPlants > 0 ) THEN
        
        IF ( time == timeNewPlants) THEN
            
            IF ( verbose ) THEN 
                WRITE(*,*) 'plants'
            END IF
        
           ! CALL PlantsGrow (t, radiation,  temperatureAir, soilMoisture, fieldCapacity, wiltingPoint,  relHumidityAir, co2 % npd (1) % value)
        
            timeNewPlants = timeNewPlants + dtPlants
            
            timeSpatialAveragePlants = time
            CALL ComputeSpatialAveragePlants  (dtPlants, lai, gpp, npp, carbonstem, carbonroot, carbonleaf, fvcover, dbh, plantsHeight, density, stemyield)
            
        END IF
        
    END IF 
    

	!-----------------------------------------------------------------------------
	!         soil water balance
	!-----------------------------------------------------------------------------
		
	IF ( dtSoilBalance > 0 ) THEN	
		IF ( time == timeNewSoilBalance ) THEN
            IF ( verbose ) THEN 
			   WRITE(*,*)'soil water balance'
            END IF
    
            CALL SolveSoilBalance ( time, raineff, irrigationFlux, &
                                    flowDirection, fvcover, vadoseDepth)
			
            timeSpatialAverageBalance = time
            CALL ComputeSpatialAverageBalance (dtSoilBalance, soilMoisture, &
                                               soilMoistureRZ, soilMoistureTZ, &
                                              runoff, infilt, percolation, &
                                              et, pet, caprise, balanceError)
				
			timeNewSoilBalance = timeNewSoilBalance + dtSoilBalance
		END IF
    ELSE
        IF ( ALLOCATED (runoff % mat) ) THEN
		  DO  i = 1, raineff % idim
			 DO j = 1, raineff % jdim
				runoff % mat(i,j) = raineff % mat(i,j)
				!percolation%mat(i,j)=0.0
			 END DO
		  END DO
        END IF
	END IF
	
	!--------------------------------------------------------------------------
	!                 groundwater
	!--------------------------------------------------------------------------
	IF (time == timeNewGroundwater) THEN
        IF ( verbose ) THEN     
	        WRITE(*,*) 'groundwater' 
        END IF
     
       
        !update river groundwater interaction discharge
        IF (riverGroundwaterInteract ) THEN
           CALL GroundwaterRiverUpdate ( waterDepth, topWidth )
        END IF

             
        !update groundwater head
        CALL GroundwaterUpdate ( time, QinSoilSub, percolation, caprise )
             
        IF ( time == timeOutGroundwater) THEN
            CALL GroundwaterOutput (time)
            timeOutGroundwater = timeOutGroundwater + dtOutGroundwater
        END IF

	    timeNewGroundwater = timeNewGroundwater + dtGroundwater
	ENDIF
       
	!-----------------------------------------------------------------------------
	!                           discharge routing
	!-----------------------------------------------------------------------------
	IF ( time == timeNewDischargeRouting ) THEN
        IF ( verbose ) THEN 
		  WRITE(*,*)'discharge routing'
        END IF
		     
        CALL DischargeRoute (dt = dtDischargeRouting, time = time, &
                        flowdir = flowDirection, &
                        runoff = runoff, dtrk = dtReservoir, &
                        riverToGroundwater = riverToGroundwater, &
                        groundwaterToriver = groundwaterToRiver, &
                        storage = storageChannelSurf)
            
 
		timeNewDischargeRouting = timeNewDischargeRouting + dtDischargeRouting   
		
  END IF
	
 
	
	!-----------------------------------------------------------------------------
	! sediment detachment 
	!-----------------------------------------------------------------------------
	
	IF (time == timeUpdateSediment) THEN
		CALL Catch ('info', 'FEST-EWB', 'update sediment erosion and deposition')
		CALL InterrillDetachmentRate (raineff)
	    timeUpdateSediment = timeUpdateSediment + dtCalcSediment
      timeSpatialAverageSediment = time
      CALL ComputeSpatialAverageSediment (dtCalcSediment, interrillErosion)
    END IF
    
    !sediment routing
    IF (routeSediment) THEN
	   IF (time == timeUpdateSedimentRouting) THEN
	     CALL SedimentRouting (dtCalcSedimentRouting, Qin, Velocita, waterDepth)


	     !update timeUpdateSedimentRouting
	     timeUpdateSedimentRouting = timeUpdateSedimentRouting + dtCalcSedimentRouting
	   END IF
	END IF

	
	!-----------------------------------------------------------------------------
	! write spatial average results
	!-----------------------------------------------------------------------------
		!output meteo spatial average
		IF (time == timeOutMeteo) THEN
            IF ( verbose ) THEN 
                WRITE(*,*) 'export meteo'
            END IF
			timeOutMeteo = timeOutMeteo + dtOutMeteo
            CALL ExportSpatialAverageMeteo ()
        END IF
        
        !output canopy interception spatial average
        IF ( time == timeOutCanopy ) THEN
            IF ( verbose ) THEN 
                WRITE(*,*) 'export canopy interception'
            END IF
            timeOutCanopy = timeOutCanopy + dtOutCanopy
            CALL ExportSpatialAverageCanopy ()
        END IF
        
        !output plants spatial average
        IF ( time == timeOutPlants ) THEN
            IF ( verbose ) THEN 
                WRITE(*,*) 'export plants'
            END IF
            timeOutPlants = timeOutPlants + dtOutPlants
            CALL ExportSpatialAveragePlants ()
        END IF
        
		!output snow spatial average
		IF (time==timeOutSnow) THEN
            IF ( verbose ) THEN 
                WRITE(*,*) 'export snow'
            END IF
			timeOutSnow = timeOutSnow + dtOutSnow
            CALL ExportSpatialAverageSnow ()
		END IF

		!output glaciers spatial average
		IF (time == timeOutIce) THEN
            IF ( verbose ) THEN 
                WRITE(*,*) 'export glaciers'
            END IF
			timeOutIce = timeOutIce + dtOutIce
            CALL ExportSpatialAverageGlaciers ()
		END IF

		!output soil balance spatial average
		IF ( time == timeOutSoilBalance ) THEN
            IF ( verbose ) THEN 
                WRITE(*,*) 'export soil balance'
            END IF
			
            CALL ExportSpatialAverageBalance ()
            
            timeOutSoilBalance = timeOutSoilBalance + dtOutSoilBalance
        END IF
		
		!!output ragguagli falda
		!IF ( time == timeOutGroundwater ) THEN
  !          IF ( verbose ) THEN 
  !              WRITE(*,*) 'export groundwater'
  !          END IF
		!	timeOutGroundwater = timeOutGroundwater+dtOutGroundwater
		!	!CALL esporta_ragguagli_falda(path_out)
		!END IF
        
		
		!output sediment spatial average
		IF (time == timeOutSediment) THEN
            IF ( verbose ) THEN 
                WRITE(*,*) 'export sediment'
            END IF
			timeOutSediment = timeOutSediment + dtOutSediment
			
        CALL ExportSpatialAverageSediment ()
			
        END IF
		!output sediment routing
		!IF (time == timeOutSedimentRouting) THEN
		!  xsOutSedimentRouting % istante = time
		!  CALL celle2stazioni(QoutSed,xsOutSedimentRouting)
		!  CALL scrivi_dato(xsOutSedimentRouting,fileunit_out_sediment_routing)
		!  !t
		!  timeOutSedimentRouting = timeOutSedimentRouting + xsOutSedimentRouting % cadenza_misure
		!
		!END IF
    
        
    !-----------------------------------------------------------------------------
	! basin average
	!-----------------------------------------------------------------------------
        
    IF ( time == timeNewBasinAverage ) THEN
        
        IF ( verbose ) THEN
            WRITE (*,*) 'export basin average'
        END IF
        
        CALL ExportBasinAverage (time, temperatureAir, &
                        temperatureAirDailyMean, temperatureAirDailyMax, &
                        temperatureAirDailyMin, precipitationRate, &
                        relHumidityAir, radiation,netRadiation, &
                        windSpeed, precipitationRateDaily, irrigationFlux, &
                        swe, soilMoisture, runoff, infilt, percolation, &
                        et, pet, deltaSoilMoisture, snowMelt, icewe, iceMelt ) 
        
        
        timeNewBasinAverage = timeNewBasinAverage + dtBasinAverage
    END IF
    
    !-----------------------------------------------------------------------------
	! raster export
	!-----------------------------------------------------------------------------
    IF ( time == timeNewRasterExport ) THEN
        
        IF ( verbose ) THEN
            WRITE (*,*) 'update raster maps'
        END IF
        
        CALL ExportMaps (time, dtCalc, temperatureAir, precipitationRate, &
                        relHumidityAir, radiation, netRadiation, &
                        windSpeed, swe, soilMoisture, runoff,  &
                        infilt, percolation, et, pet ) 
        
        
        timeNewRasterExport = timeNewRasterExport + dtCalc
    END IF
    
    
    
        
    !-----------------------------------------------------------------------------
	! save status
	!-----------------------------------------------------------------------------
    IF (  hotstart ) THEN
        !IF ( IsTimeHotStart (time) ) THEN
        IF ( CronIsTime (time, cronHotstart) ) THEN
            IF ( verbose ) THEN
		         WRITE (*,*) 'save state variables for hot start'
            END IF
		    CALL SaveHotStart (time)
        END IF
    END IF
    
	!-----------------------------------------------------------------------------
	! time forward
	!-----------------------------------------------------------------------------
	time = time + dtCalc
	
	
	IF ( time > timeStop ) THEN
		IF (  hotstart ) THEN
            IF ( saveLast > 0 ) THEN
                IF ( verbose ) THEN
		           WRITE (*,*) 'save final state variables for hot start'
                END IF
			    CALL SaveHotStart
            END IF
        END IF
	END IF
END DO

!final message
IF ( verbose ) THEN 
    WRITE(*,*) 'simulation finished'
END IF

!-----------------------------------------
!    terminate log
!-----------------------------------------

CALL LogStop ()


CONTAINS
    
!==============================================================================
!| Description:
!  Save state variables for starting a new simulation from the current condition
!  
SUBROUTINE SaveHotStart &
!
(time)

IMPLICIT NONE


TYPE(DateTime), OPTIONAL :: time
CHARACTER (LEN = 26)   :: prefix
CHARACTER (LEN = 300)  :: file


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

IF (PRESENT(time)) THEN
	prefix = time
	prefix = prefix (1:19) // '_'
	prefix (14:14) = '-'
	prefix (17:17) = '-'
		
ELSE
	prefix = '                    '
END IF

!snow
IF ( dtSnow > 0 ) THEN
    
    !swe
    file = pathHotstart (1:LEN_TRIM (pathHotstart)) // &
	                  prefix (1:LEN_TRIM(prefix)) // 'swe'
    
    CALL ExportGrid (layer = swe, fileName = file, &
                     fileFormat = ESRI_BINARY)
    
    !water in snow
    file = pathHotstart (1:LEN_TRIM (pathHotstart)) // &
	                  prefix (1:LEN_TRIM(prefix)) // 'water-snow'
    
    CALL ExportGrid (layer = waterInSnow, fileName = file, &
                     fileFormat = ESRI_BINARY)

END IF

!glaciers
IF ( dtIce > 0 ) THEN
    
    !ice water equivalent
    file = pathHotstart (1:LEN_TRIM (pathHotstart)) // &
	                  prefix (1:LEN_TRIM(prefix)) // 'icewe'
    
    CALL ExportGrid (layer = icewe, fileName = file, &
                     fileFormat = ESRI_BINARY)
    
    !water in ice
    file = pathHotstart (1:LEN_TRIM (pathHotstart)) // &
	                  prefix (1:LEN_TRIM(prefix)) // 'water-ice'
    
    CALL ExportGrid (layer = waterInIce, fileName = file, &
                     fileFormat = ESRI_BINARY)

END IF

!soil balance
IF ( dtSoilBalance > 0 ) THEN

    !soil saturation root zone
    file = pathHotstart (1:LEN_TRIM (pathHotstart)) // &
	                  prefix (1:LEN_TRIM(prefix)) // 'sat-rz'
    
    CALL ExportGrid (layer = soilSatRZ, fileName = file, &
                     fileFormat = ESRI_BINARY)
    
    !soil saturation transmission zone
    file = pathHotstart (1:LEN_TRIM (pathHotstart)) // &
	                  prefix (1:LEN_TRIM(prefix)) // 'sat-tz'
    
    CALL ExportGrid (layer = soilSatTZ, fileName = file, &
                    fileFormat = ESRI_BINARY)
    
    !precipitation status
    file = pathHotstart (1:LEN_TRIM (pathHotstart)) // &
	                  prefix (1:LEN_TRIM(prefix)) // 'rainflag'
    
    CALL ExportGrid (layer = rainFlag, fileName = file, &
                    fileFormat = ESRI_BINARY)
    
    !interstorm duration
    file = pathHotstart (1:LEN_TRIM (pathHotstart)) // &
	                  prefix (1:LEN_TRIM(prefix)) // 'interstorm'
    
    CALL ExportGrid (layer = interstormDuration, fileName = file, &
                    fileFormat = ESRI_BINARY)
    
    IF ( infiltrationModel == SCS_CN ) THEN
        
        !soil retention
        file = pathHotstart (1:LEN_TRIM (pathHotstart)) // &
	                  prefix (1:LEN_TRIM(prefix)) // 'seff'
    
        CALL ExportGrid (layer = sEff, fileName = file, &
                         fileFormat = ESRI_BINARY)
        
        !cumulative precipitation
        file = pathHotstart (1:LEN_TRIM (pathHotstart)) // &
	                  prefix (1:LEN_TRIM(prefix)) // 'raincum'
    
        CALL ExportGrid (layer = raincum, fileName = file, &
                         fileFormat = ESRI_BINARY)
    END IF
    
    IF ( infiltrationModel == PHILIPEQ .OR. &
         infiltrationModel == GREEN_AMPT ) THEN
        
        !cumulative infiltration
        file = pathHotstart (1:LEN_TRIM (pathHotstart)) // &
	                  prefix (1:LEN_TRIM(prefix)) // 'cuminf'
    
        CALL ExportGrid (layer = cuminf, fileName = file, &
                         fileFormat = ESRI_BINARY)
    END IF
END IF

! discharge routing
IF ( dtDischargeRouting > 0 ) THEN
    !Qin
    file = pathHotstart (1:LEN_TRIM (pathHotstart)) // &
	                  prefix (1:LEN_TRIM(prefix)) // 'Qin'
    
    CALL ExportGrid (layer = Qin, fileName = file, &
                     fileFormat = ESRI_BINARY)
    
    !Qout
    file = pathHotstart (1:LEN_TRIM (pathHotstart)) // &
	                  prefix (1:LEN_TRIM(prefix)) // 'Qout'
    
    CALL ExportGrid (layer = Qout, fileName = file, &
                     fileFormat = ESRI_BINARY)
    
    !Qlat
    file = pathHotstart (1:LEN_TRIM (pathHotstart)) // &
	                  prefix (1:LEN_TRIM(prefix)) // 'Qlat'
    
    CALL ExportGrid (layer = Plat, fileName = file, &
                     fileFormat = ESRI_BINARY)
    ! diversions
    IF ( dtDiversion > 0 ) THEN
        IF ( PRESENT (time) ) THEN
           CALL DiversionSaveStatus ( pathHotstart, time )
        ELSE
           CALL DiversionSaveStatus ( pathHotstart )
        END IF
    END IF
    
    ! reservoirs
    IF ( dtReservoir > 0 ) THEN
        IF ( PRESENT (time) ) THEN
           CALL ReservoirSaveStatus ( pathHotstart, time )
        ELSE
           CALL ReservoirSaveStatus ( pathHotstart )
        END IF
    END IF
END IF

!	!Falda
!IF (dtGroundwater > 0) THEN
!	DO I=1,AltezzaPiezometrica%NumeroStrati             
!		                           
!	     CALL write_raster_bin(path_hotstart(1:LEN_TRIM(path_hotstart))// &
!	                  pref_s(1:LEN_TRIM(pref_s))//'strato'//CHAR(64+I), &
!	                  matrice_real=AltezzaPiezometrica%piezometria(I))
!	                                           
!		 IF (interazione_falda_fiume ) THEN
!		  
!		    CALL write_raster_bin(path_hotstart(1:LEN_TRIM(path_hotstart))// &
!	                  pref_s(1:LEN_TRIM(pref_s))//'baseflow_falda', &
!	                  matrice_real=BaseFlow_falda)
!	                  
!	        CALL write_raster_bin(path_hotstart(1:LEN_TRIM(path_hotstart))// &
!	                  pref_s(1:LEN_TRIM(pref_s))//'ricarica_falda', &
!	                  matrice_real=ricarica_falda_fiume)
!	                                            
!		 END IF
!	ENDDO
!END IF
	
!sediment
!IF (dtCalcSediment > 0) THEN
!  CALL write_raster_bin(path_hotstart(1:LEN_TRIM(path_hotstart))// &
!	                  pref_s(1:LEN_TRIM(pref_s))//'deltaSed',matrice_real=deltaSed)
!END IF
!
!IF (dtCalcSedimentRouting > 0) THEN
!
!!to do
!
!END IF
	


END SUBROUTINE SaveHotStart



SUBROUTINE Splash ()

WRITE (*,*) 


!WRITE (*,'(a)') '  ______ ______  _____ _______ '
!WRITE (*,'(a)') ' |  ____|  ____|/ ____|__   __|'
!WRITE (*,'(a)') ' | |__  | |__  | (___    | |   '
!WRITE (*,'(a)') ' |  __| |  __|  \___ \   | |   '
!WRITE (*,'(a)') ' | |    | |____ ____) |  | |   '
!WRITE (*,'(a)') ' |_|    |______|_____/   |_|   '

!WRITE (*,'(a)') 
!WRITE (*,'(a)') '-o--o o--o  o-o  o--O-o '
!WRITE (*,'(a)') ' |    |    |        |   '
!WRITE (*,'(a)') ' O-o  O-o   o-o     |   '
!WRITE (*,'(a)') ' |    |        |    |   '
!WRITE (*,'(a)') ' o    o--o o--o     o   '
!WRITE (*,'(a)')      

                
WRITE (*,'(a)') ' ________  _______       ________       _________    '
WRITE (*,'(a)') '|\  _____\|\  ___ \     |\   ____\     |\___   ___\  '
WRITE (*,'(a)') '\ \  \__/ \ \   __/|    \ \  \___|_    \|___ \  \_|  '
WRITE (*,'(a)') ' \ \   __\ \ \  \_|/__   \ \_____  \        \ \  \   '
WRITE (*,'(a)') '  \ \  \_|  \ \  \_|\ \   \|____|\  \        \ \  \  '
WRITE (*,'(a)') '   \ \__\    \ \_______\    ____\_\  \        \ \__\ '
WRITE (*,'(a)') '    \|__|     \|_______|   |\_________\        \|__| '
WRITE (*,'(a)') '                           \|_________|              '
WRITE (*,'(a)') '                                                     '
! (*,'(a)') '        E N E R G Y - W A T E R - B A L A N C E      '
WRITE (*,'(a)') '     S P A T I A L L Y - D I S T R I B U T E D '
WRITE (*,'(a)') '        H Y D R O L O G I C A L - M O D E L    '


!WRITE (*,'(a)') ' Flood Event Spatially distributed rainfall-runoff Transformation model'
WRITE (*,'(a)') 

RETURN
END SUBROUTINE Splash






!-------------------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------------------
END PROGRAM Fest