!! Compute average value data within a river basin 
!|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.2 - 25th October 2024    
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 19/Jan/2024 | Original code |
! | 1.1      | 25/Mar/2024 | added swe and soil balance terms |
! | 1.2      | 25/Oct/2024 | added snow and ice melt and ice water equivalent |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! routines to compute average value of grid data over a given river basin 
! Firstly variable is cumulated along the river network 
! by following the flow direction. Then the cumulated value
! in a given point divided by the number of cell from
! the flow accumulation map is written to output file. 
! The user provides a list of point coordinates 
! that define the closing section of the river basins
! and set the variables to be processed in the configuration file
! like in the following example
!```
! # set the variables to export
! # 1 activates export
! # 0 suppresses export
!
! # meteo
! precipitation = 1
! daily-precipitation = 0
! temperature = 1
! temperature-daily-mean = 0
! temperature-daily-max = 0
! temperature-daily-min = 0
! relative-humidity = 0
! solar-radiation = 0
! net-radiation = 0
! wind-speed = 0
! irrigation = 0
!
! #snow
! snow-water-equivalent = 1
! snow-melt = 1
!
! #glaciers
! ice-water-equivalent = 0
! ice-melt = 0
!
! #soil balance
! soil-moisture = 1
! delta-soil-moisture = 1
! runoff = 1
! infiltration = 1
! percolation = 1
! actual-ET = 1
! potential-ET = 0
!
!```
!   
! The value is computed for all variables marked by 1. 
! When one variable is marked by 1 but it is not allocated 
! because not computed by the FEST model according to options
! defined in the configuration files, value is not exported. 
! For example, if user set to export wind-speed 
! but windspeed is not used in
! the current simulation, values of windspeed are 
! not written in the output file.
! One output file is created for each variable.
! So average temperature values for all points 
! are written in a file, precipitation values are written in 
! a different file, and so on.
!
! Variables that can be exported, description, and unit are listed 
! in the following table.
!
! | variable               | Description                                                        | Unit           |
! |------------------------|--------------------------------------------------------------------|----------------|
! | precipitation          | Precipitation fallen in the current time step                      | mm             |
! | daily-precipitation    | Precipitation fallen in 24 hours                                   | mm             |
! | temperature            | Air temperature of the current time step fallen in 24 hours        | Celsius degree |
! | temperature-daily-mean | Mean daily air temperature                                         | Celsius degree |
! | temperature-daily-max  | Maximum daily air temperature                                      | Celsius degree |
! | temperature-daily-min  | Maximum daily air temperature                                      | Celsius degree |
! | relative-humidity      | Air relative humidity of the current time step                     | % (0-100)      |
! | solar-radiation        | Solar radiation of the current time step                           | w/m²           |
! | net-radiation          | Net radiation of the current time step                             | w/m²           |
! | wind-speed             | Wind speed of the current time step                                | m/s            |
! | irrigation             | Irrigation amount of the current time step                         | mm             |
! | snow-water-equivalent  | Snow water equivalent of the current time step                     | mm             |
! | snow-melt              | Snow melt of the current time step                                 | mm             |
! | ice-water-equivalent   | Glaciers ice water equivalent of the current time step             | mm             |
! | ice-melt               | Glaciers melt of the current time step                             | mm             |
! | soil-moisture          | Soil moisture of the current time step                             | \-             |
! | runoff                 | Runoff of the current time step                                    | mm             |
! | infiltration           | Infiltration into soil of the current time step                    | mm             |
! | percolation            | Deep percolation out of transmission zone of the current time step | mm             |
! | actual-ET              | Actual evapotranspiration of the current time step                 | mm             |
! | potential-ET           | Potential evapotranspiration of the current time step              | mm             |
! | delta-soil-moisture    | Change in soil water storage of the current time step              | mm             |
!
!
! The name of output files is the concatenation of result 
! folder name defined in the main configuration file <folder>, 
! and a suffix that reminds the name of variable, 
! as listed in the following table.
!
! | variable                    | Output file name                          |
! |-----------------------------|-------------------------------------------|
! | precipitation               | `<folder>` `mean_precipitation.fts`       |
! | daily-precipitation         | `<folder>` `mean_pdaily.fts`              |
! | temperature                 | `<folder>` `mean_temperature.fts`         |
! | temperature-daily-mean      | `<folder>` `mean_tmean.fts`               |
! | temperature-daily-max       | `<folder>` `mean_tmax.fts`                |
! | temperature-daily-min       | `<folder>` `mean_tmin.fts`                |
! | relative-humidity           | `<folder>` `mean_rh.fts`                  |
! | solar-radiation             | `<folder>` `mean_rad.fts`                 |
! | net-radiation               | `<folder>` `mean_netrad.fts`              |
! | wind-speed                  | `<folder>` `mean_windspeed.fts`           |
! | irrigation                  | `<folder>` `mean_irrigation.fts`          |
! | snow-water-equivalent       | `<folder>` `mean_swe.fts`                 |
! | snow-melt                   | `<folder>` `mean_snow-melt.fts`           |
! | ice-water-equivalent        | `<folder>` `mean_icewe.fts`               |
! | ice-melt                    | `<folder>` `mean_ice-melt.fts`            |
! | soil-moisture               | `<folder>` `mean_soil-moisture.fts`       |
! | runoff                      | `<folder>` `mean_runoff.fts`              |
! | infiltration                | `<folder>` `mean_infiltration.fts`        |
! | percolation                 | `<folder>` `mean_percolation.fts`         |
! | actual-ET                   | `<folder>` `mean_et.fts`                  |
! | potential-ET                | `<folder>` `mean_pet.fts`                 |
! | delta-soil-moisture         | `<folder>` `mean_delta-soil-moisture.fts` |
!
MODULE BasinAverage

! Modules used: 

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

USE DomainProperties, ONLY : &
    !imported variables:
    mask

USE MorphologicalProperties, ONLY : &
    !Imported variables:
    flowDirection, &
    flowAccumulation, &
    streamNetwork, &
    flowDirection_loaded, &
    flowAccumulation_loaded, &
    streamNetworkCreated

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

USE ObservationalNetworks, ONLY : &
    !imported definitions:
    ObservationalNetwork, &
    !Imported routines:
    ReadMetadata, &
    CopyNetwork, &
    WriteMetaData, &
    WriteData, &
    AssignDataFromGrid

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

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

USE Chronos, ONLY : &
    !Imported definitions:
    DateTime, &
    !Imported operands:
    ASSIGNMENT( = )

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

USE Loglib, ONLY : &
    !Imported routines:
    Catch

USE Utilities, ONLY : &
    !imported routines:
    GetUnit

IMPLICIT NONE 

INTEGER (KIND = short) :: dtBasinAverage

!Public routines
PUBLIC :: InitBasinAverage
PUBLIC :: ReadPointFileBasinAverage
PUBLIC :: ExportBasinAverage


!private declarations
INTEGER (KIND = short), PRIVATE :: countVar  !!count number of variables active for output 
TYPE (ObservationalNetwork), PRIVATE :: sites
TYPE (ObservationalNetwork), PRIVATE :: sitesPrecipitation
TYPE (ObservationalNetwork), PRIVATE :: sitesPrecipitationDaily
TYPE (ObservationalNetwork), PRIVATE :: sitesTemp
TYPE (ObservationalNetwork), PRIVATE :: sitesTmean
TYPE (ObservationalNetwork), PRIVATE :: sitesTmin
TYPE (ObservationalNetwork), PRIVATE :: sitesTmax
TYPE (ObservationalNetwork), PRIVATE :: sitesRH
TYPE (ObservationalNetwork), PRIVATE :: sitesRadiation
TYPE (ObservationalNetwork), PRIVATE :: sitesNetRadiation
TYPE (ObservationalNetwork), PRIVATE :: sitesWindSpeed
TYPE (ObservationalNetwork), PRIVATE :: sitesIrrigation
TYPE (ObservationalNetwork), PRIVATE :: sitesSWE
TYPE (ObservationalNetwork), PRIVATE :: sitesSnowMelt
TYPE (ObservationalNetwork), PRIVATE :: sitesIWE
TYPE (ObservationalNetwork), PRIVATE :: sitesIceMelt
TYPE (ObservationalNetwork), PRIVATE :: sitesSM
TYPE (ObservationalNetwork), PRIVATE :: sitesDSM
TYPE (ObservationalNetwork), PRIVATE :: sitesRunoff
TYPE (ObservationalNetwork), PRIVATE :: sitesInfiltration
TYPE (ObservationalNetwork), PRIVATE :: sitesPercolation
TYPE (ObservationalNetwork), PRIVATE :: sitesET
TYPE (ObservationalNetwork), PRIVATE :: sitesPET

INTEGER (KIND = short), PRIVATE :: fileUnitPrecipitation
INTEGER (KIND = short), PRIVATE :: fileUnitPrecipitationDaily
INTEGER (KIND = short), PRIVATE :: fileUnitTemp
INTEGER (KIND = short), PRIVATE :: fileUnitTmean
INTEGER (KIND = short), PRIVATE :: fileUnitTmin
INTEGER (KIND = short), PRIVATE :: fileUnitTmax
INTEGER (KIND = short), PRIVATE :: fileUnitRH
INTEGER (KIND = short), PRIVATE :: fileUnitRadiation
INTEGER (KIND = short), PRIVATE :: fileUnitNetRadiation
INTEGER (KIND = short), PRIVATE :: fileUnitWindSpeed
INTEGER (KIND = short), PRIVATE :: fileUnitIrrigation
INTEGER (KIND = short), PRIVATE :: fileUnitSWE
INTEGER (KIND = short), PRIVATE :: fileUnitSnowMelt
INTEGER (KIND = short), PRIVATE :: fileUnitIWE
INTEGER (KIND = short), PRIVATE :: fileUnitIceMelt
INTEGER (KIND = short), PRIVATE :: fileUnitSM
INTEGER (KIND = short), PRIVATE :: fileUnitDSM
INTEGER (KIND = short), PRIVATE :: fileUnitRunoff
INTEGER (KIND = short), PRIVATE :: fileUnitInfiltration
INTEGER (KIND = short), PRIVATE :: fileUnitPercolation
INTEGER (KIND = short), PRIVATE :: fileUnitET
INTEGER (KIND = short), PRIVATE :: fileUnitPET

!temporary grid to cumulate data
TYPE (grid_real), PRIVATE :: cum


!active output
LOGICAL, PRIVATE :: varOut (22)  !1 = precipitation, 
                        !2 = daily-precipitation,
                        !3 = air-temperature, 
                        !4 = air-temperature-daily-mean,
                        !5 = air-temperature-daily-max
                        !6 = air-temperature-daily-min
                        !7 = relative-humidity
                        !8 = solar-radiation,
                        !9 = net-radiation
                        !10 = wind-speed
                        !11 = irrigation
                        !12 = snow-water-equivalent
                        !13 = soil-moisture
                        !14 = runoff
                        !15 = infiltration
                        !16 = percolation
                        !17 = actual-ET
                        !18 = potential-ET
                        !19 = delta-soil-moisture
                        !20 = snow-melt
                        !21 = ice-water-equivalent
                        !22 = ice-melt


!private routines
PRIVATE :: BasinCumulate
PRIVATE :: BasinRasterExport



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


!==============================================================================
!| Description: 
!  Read point file
SUBROUTINE ReadPointFileBasinAverage   & 
!
 ( string )  

IMPLICIT NONE

!arguments with intent in:
CHARACTER(LEN = *), INTENT(IN)    :: string 

!local declarations

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

CALL ReadMetadata ( string, sites )

dtBasinAverage = sites % timeIncrement

RETURN
END SUBROUTINE ReadPointFileBasinAverage
 

    

!==============================================================================
!| Description: 
!  Initialization of basin average
SUBROUTINE InitBasinAverage   & 
!
 (fileini, pathout, temp, tmean, tmax, tmin, precipitation, &
  rh, radiation, netradiation, windspeed, daily_precipitation, &
  irrigation, swe, sm, runoff, infiltration, percolation, et, &
  pet, dsm, snowmelt, iwe, icemelt  )  

IMPLICIT NONE

!arguments with intent in:
CHARACTER(LEN = *), INTENT(IN)    :: fileini 
CHARACTER(LEN = *), INTENT(IN)    :: pathout   
TYPE (grid_real), INTENT(IN) :: temp !!air temperarure (°C)
TYPE (grid_real), INTENT(IN) :: tmean !!air temperarure daily mean(°C)
TYPE (grid_real), INTENT(IN) :: tmax !!air temperarure daily max (°C)
TYPE (grid_real), INTENT(IN) :: tmin !!air temperarure daily min (°C)
TYPE (grid_real), INTENT(IN) :: precipitation !!precipitation rate (m/s)
TYPE (grid_real), INTENT(IN) :: rh !!air relative humidity (0-100)
TYPE (grid_real), INTENT(IN) :: radiation !!solar radiation (w/m2)
TYPE (grid_real), INTENT(IN) :: netradiation !!net radiation (w/m2)
TYPE (grid_real), INTENT(IN) :: windspeed !!wind speed (m/s)
TYPE (grid_real), INTENT(IN) :: daily_precipitation !!daily precipitation rate (m/s)
TYPE (grid_real), INTENT(IN) :: irrigation !!irrigation rate (m/s)
TYPE (grid_real), INTENT(IN) :: swe !!snow water equivalent (m)
TYPE (grid_real), INTENT(IN) :: sm !!soil mositure (-)
TYPE (grid_real), INTENT(IN) :: runoff !!runoff (m/s)
TYPE (grid_real), INTENT(IN) :: infiltration !!infiltration (m/s)
TYPE (grid_real), INTENT(IN) :: percolation !!deep percolation (m/s)
TYPE (grid_real), INTENT(IN) :: et !!actual evapotranspiration (m/s)
TYPE (grid_real), INTENT(IN) :: pet !!potential evapotranspiration (m/s)
TYPE (grid_real), INTENT(IN) :: dsm !!soil moisture variation (m)
TYPE (grid_real), INTENT(IN) :: snowmelt !! snow melt (m)
TYPE (grid_real), INTENT(IN) :: iwe !! ice water equivalent (m)
TYPE (grid_real), INTENT(IN) :: icemelt !! ice melt (m)

 

!local declarations
TYPE(IniList)          :: iniDB
CHARACTER (LEN = 300)  :: string
!-------------------------------end of declaration-----------------------------


!Check morphological properties are available

IF ( .NOT. flowDirection_loaded ) THEN
    CALL Catch ('error', 'BasinAverage', 'flow direction missing, &
                                  check morphological properties')
END IF

IF ( .NOT. flowAccumulation_loaded ) THEN
    CALL Catch ('error', 'BasinAverage', 'flow accumulation missing, &
                                  check morphological properties')
END IF

IF ( .NOT. streamNetworkCreated ) THEN
    CALL Catch ('error', 'BasinAverage', 'stream network missing, &
                                  check morphological properties')
END IF

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

!create and export basin raster maps
IF (KeyIsPresent ('raster-export', iniDB) ) THEN
    IF ( IniReadInt ('raster-export', iniDB)  == 1 ) THEN
        CALL BasinRasterExport ( pathout )
    END IF
END IF


! search for active variable for output
CALL Catch ('info', 'BasinAverage', 'checking for active variables ')

countVar = 0
!precipitation
IF ( IniReadInt ('precipitation', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (precipitation % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', 'precipitation not allocated, &
                                        forced to not export basin average ')
       varOut (1) = .FALSE.
   ELSE
       varOut (1) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesPrecipitation )
       fileUnitPrecipitation = GetUnit ()
       string = TRIM(pathout) // 'mean_precipitation.fts'
       OPEN ( unit = fileUnitPrecipitation, file = string )
       sitesPrecipitation % description = 'average precipitation'
       sitesPrecipitation % unit = 'mm'
       sitesPrecipitation % offsetZ = 0.
       CALL WriteMetadata ( network = sitesPrecipitation, &
                        fileunit = fileUnitPrecipitation )
       CALL WriteData (sitesPrecipitation, fileUnitPrecipitation, .TRUE.) 
       
   END IF
ELSE
   varOut (1) = .FALSE.
END IF

!daily-precipitation
IF ( IniReadInt ('daily-precipitation', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (daily_precipitation % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', &
                   'daily_precipitation not allocated, &
                    forced to not export basin average ')
       varOut (2) = .FALSE.
   ELSE
       varOut (2) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesPrecipitationDaily )
       fileUnitPrecipitationDaily = GetUnit ()
       string = TRIM(pathout) // 'mean_pdaily.fts'
       OPEN ( unit = fileUnitPrecipitationDaily, file = string )
       sitesPrecipitationDaily % description = 'daily precipitation'
       sitesPrecipitationDaily % unit = 'mm'
       sitesPrecipitationDaily % offsetZ = 0.
       CALL WriteMetadata ( network = sitesPrecipitationDaily, &
                        fileunit = fileUnitPrecipitationDaily )
       CALL WriteData (sitesPrecipitationDaily, &
                       fileUnitPrecipitationDaily, .TRUE.) 
   END IF
ELSE
   varOut (2) = .FALSE.
END IF

!air-temperature
IF ( IniReadInt ('temperature', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (temp % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', &
                   'air temperature not allocated, &
                   forced to not export basin average ')
       varOut (3) = .FALSE.
   ELSE
       varOut (3) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesTemp )
       fileUnitTemp = GetUnit ()
       string = TRIM(pathout) // 'mean_temperature.fts'
       OPEN ( unit = fileUnitTemp, file = string )
       sitesTemp % description = 'air temperature'
       sitesTemp % unit = 'Degree Celsius'
       sitesTemp % offsetZ = 0.
       CALL WriteMetadata ( network = sitesTemp, &
                        fileunit = fileUnitTemp )
       CALL WriteData (sitesTemp, fileUnitTemp, .TRUE.)
   END IF
ELSE
   varOut (3) = .FALSE.
END IF

!temperature-daily-mean
IF ( IniReadInt ('temperature-daily-mean', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (tmean % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', &
                   'mean daily air temperature not allocated, &
                    forced to not export basin average ')
       varOut (4) = .FALSE.
   ELSE
       varOut (4) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesTmean )
       fileUnitTmean = GetUnit ()
       string = TRIM(pathout) // 'mean_tmean.fts'
       OPEN ( unit = fileUnitTmean, file = string )
       sitesTmean % description = 'mean daily air temperature'
       sitesTmean % unit = 'Degree Celsius'
       sitesTmean % offsetZ = 0.
       CALL WriteMetadata ( network = sitesTmean, &
                        fileunit = fileUnitTmean )
       CALL WriteData (sitesTmean, fileUnitTmean, .TRUE.)
   END IF
ELSE
   varOut (4) = .FALSE.
END IF

!air-temperature-daily-max
IF ( IniReadInt ('temperature-daily-max', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (tmax % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', &
                   'max daily air temperature not allocated, &
                    forced to not export basin average ')
       varOut (5) = .FALSE.
   ELSE
       varOut (5) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesTmax )
       fileUnitTmax = GetUnit ()
       string = TRIM(pathout) // 'mean_tmax.fts'
       OPEN ( unit = fileUnitTmax, file = string )
       sitesTmean % description = 'maximum daily air temperature'
       sitesTmean % unit = 'Degree Celsius'
       sitesTmean % offsetZ = 0.
       CALL WriteMetadata ( network = sitesTmax, &
                        fileunit = fileUnitTmax )
       CALL WriteData (sitesTmax, fileUnitTmax, .TRUE.)
   END IF
ELSE
   varOut (5) = .FALSE.
END IF


!air-temperature-daily-min
IF ( IniReadInt ('temperature-daily-min', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (tmin % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', &
                   'min daily air temperature not allocated, &
                    forced to not export basin average ')
       varOut (6) = .FALSE.
   ELSE
       varOut (6) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesTmin )
       fileUnitTmin = GetUnit ()
       string = TRIM(pathout) // 'mean_tmin.fts'
       OPEN ( unit = fileUnitTmin, file = string )
       sitesTmin % description = 'minimum daily air temperature'
       sitesTmean % unit = 'Degree Celsius'
       sitesTmean % offsetZ = 0.
       CALL WriteMetadata ( network = sitesTmin, &
                        fileunit = fileUnitTmin )
       CALL WriteData (sitesTmin, fileUnitTmin, .TRUE.)
   END IF
ELSE
   varOut (6) = .FALSE.
END IF


!relative-humidity
IF ( IniReadInt ('relative-humidity', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (rh % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', &
                   'relative humidity not allocated, &
                    forced to not export basin average ')
       varOut (7) = .FALSE.
   ELSE
       varOut (7) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesRH )
       fileUnitRH = GetUnit ()
       string = TRIM(pathout) // 'mean_rh.fts'
       OPEN ( unit = fileUnitRH, file = string )
       sitesRH % description = 'mean air relative humidity'
       sitesRH % unit = '% 0-100'
       sitesRH % offsetZ = 0.
       CALL WriteMetadata ( network = sitesRH, &
                        fileunit = fileUnitRH )
       CALL WriteData (sitesRH, fileUnitRH, .TRUE.)
   END IF
ELSE
   varOut (7) = .FALSE.
END IF

!solar-radiation
IF ( IniReadInt ('solar-radiation', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (radiation % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', &
                   'solar radiation not allocated, &
                    forced to not export basin average ')
       varOut (8) = .FALSE.
   ELSE
       varOut (8) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesRadiation )
       fileUnitRadiation = GetUnit ()
       string = TRIM(pathout) // 'mean_rad.fts'
       OPEN ( unit = fileUnitRadiation, file = string )
       sitesRadiation % description = 'mean solar radiation'
       sitesRadiation % unit = 'w/m2'
       sitesRadiation % offsetZ = 0.
       CALL WriteMetadata ( network = sitesRadiation, &
                        fileunit = fileUnitRadiation )
       CALL WriteData (sitesRadiation, fileUnitRadiation, .TRUE.)
   END IF
ELSE
   varOut (8) = .FALSE.
END IF

!net-radiation
IF ( IniReadInt ('net-radiation', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (netradiation % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', 'net radiation not allocated, &
                                        forced to not export basin average ')
       varOut (9) = .FALSE.
   ELSE
       varOut (9) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesNetRadiation )
       fileUnitNetRadiation = GetUnit ()
       string = TRIM(pathout) // 'mean_netrad.fts'
       OPEN ( unit = fileUnitNetRadiation, file = string )
       sitesNetRadiation % description = 'mean net radiation'
       sitesNetRadiation % unit = 'w/m2'
       sitesNetRadiation % offsetZ = 0.
       CALL WriteMetadata ( network = sitesNetRadiation, &
                        fileunit = fileUnitNetRadiation )
       CALL WriteData (sitesNetRadiation, fileUnitNetRadiation, .TRUE.)
   END IF
ELSE
   varOut (9) = .FALSE.
END IF

!wind-speed
IF ( IniReadInt ('wind-speed', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (windspeed % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', 'wind speed not allocated, &
                                        forced to not export basin average ')
       varOut (10) = .FALSE.
   ELSE
       varOut (10) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesWindSpeed )
       fileUnitWindSpeed = GetUnit ()
       string = TRIM(pathout) // 'mean_windspeed.fts'
       OPEN ( unit = fileUnitWindSpeed, file = string )
       sitesWindSpeed % description = 'mean wind speed'
       sitesWindSpeed % unit = 'm/s'
       sitesWindSpeed % offsetZ = 0.
       CALL WriteMetadata ( network = sitesWindSpeed, &
                        fileunit = fileUnitWindSpeed )
       CALL WriteData (sitesWindSpeed, fileUnitWindSpeed, .TRUE.)
   END IF
ELSE
   varOut (10) = .FALSE.
END IF

!irrigation
IF ( IniReadInt ('irrigation', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (irrigation % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', 'irrigation not allocated, &
                                        forced to not export basin average ')
       varOut (11) = .FALSE.
   ELSE
       varOut (11) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesIrrigation )
       fileUnitIrrigation = GetUnit ()
       string = TRIM(pathout) // 'mean_irrigation.fts'
       OPEN ( unit = fileUnitIrrigation, file = string )
       sitesIrrigation % description = 'mean irrigation'
       sitesIrrigation % unit = 'mm'
       sitesIrrigation % offsetZ = 0.
       CALL WriteMetadata ( network = sitesIrrigation, &
                        fileunit = fileUnitIrrigation )
       CALL WriteData (sitesIrrigation, fileUnitIrrigation, .TRUE.)
   END IF
ELSE
   varOut (11) = .FALSE.
END IF

!snow water equivalent
IF ( IniReadInt ('snow-water-equivalent', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (swe % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', 'SWE not allocated, &
                                        forced to not export basin average ')
       varOut (12) = .FALSE.
   ELSE
       varOut (12) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesSWE )
       fileUnitSWE = GetUnit ()
       string = TRIM(pathout) // 'mean_swe.fts'
       OPEN ( unit = fileUnitSWE, file = string )
       sitesSWE % description = 'mean snow water equivalent'
       sitesSWE % unit = 'mm'
       sitesSWE % offsetZ = 0.
       CALL WriteMetadata ( network = sitesSWE, &
                        fileunit = fileUnitSWE )
       CALL WriteData (sitesSWE, fileUnitSWE, .TRUE.)
   END IF
ELSE
   varOut (12) = .FALSE.
END IF

!soil moisture
IF ( IniReadInt ('soil-moisture', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (sm % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', 'soil-moisture not allocated, &
                                        forced to not export basin average ')
       varOut (13) = .FALSE.
   ELSE
       varOut (13) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesSM )
       fileUnitSM = GetUnit ()
       string = TRIM(pathout) // 'mean_soil-moisture.fts'
       OPEN ( unit = fileUnitSM, file = string )
       sitesSM % description = 'mean soil moisture'
       sitesSM % unit = 'm3/m3'
       sitesSM % offsetZ = 0.
       CALL WriteMetadata ( network = sitesSM, &
                        fileunit = fileUnitSM )
       CALL WriteData (sitesSM, fileUnitSM, .TRUE.)
   END IF
ELSE
   varOut (13) = .FALSE.
END IF

!runoff
IF ( IniReadInt ('runoff', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (runoff % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', 'runoff not allocated, &
                                        forced to not export basin average ')
       varOut (14) = .FALSE.
   ELSE
       varOut (14) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesRunoff )
       fileUnitRunoff = GetUnit ()
       string = TRIM(pathout) // 'mean_runoff.fts'
       OPEN ( unit = fileUnitRunoff, file = string )
       sitesRunoff % description = 'mean runoff'
       sitesRunoff % unit = 'mm'
       sitesRunoff % offsetZ = 0.
       CALL WriteMetadata ( network = sitesRunoff, &
                        fileunit = fileUnitRunoff )
       CALL WriteData (sitesRunoff, fileUnitRunoff, .TRUE.)
   END IF
ELSE
   varOut (14) = .FALSE.
END IF


!infiltration
IF ( IniReadInt ('infiltration', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (infiltration % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', 'infiltration not allocated, &
                                        forced to not export basin average ')
       varOut (15) = .FALSE.
   ELSE
       varOut (15) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesInfiltration )
       fileUnitInfiltration = GetUnit ()
       string = TRIM(pathout) // 'mean_infiltration.fts'
       OPEN ( unit = fileUnitInfiltration, file = string )
       sitesInfiltration % description = 'mean infiltration'
       sitesInfiltration % unit = 'mm'
       sitesInfiltration % offsetZ = 0.
       CALL WriteMetadata ( network = sitesInfiltration, &
                        fileunit = fileUnitInfiltration )
       CALL WriteData (sitesInfiltration, fileUnitInfiltration, .TRUE.)
   END IF
ELSE
   varOut (15) = .FALSE.
END IF


!percolation
IF ( IniReadInt ('percolation', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (percolation % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', 'percolation not allocated, &
                                        forced to not export basin average ')
       varOut (16) = .FALSE.
   ELSE
       varOut (16) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesPercolation )
       fileUnitPercolation = GetUnit ()
       string = TRIM(pathout) // 'mean_percolation.fts'
       OPEN ( unit = fileUnitPercolation, file = string )
       sitesPercolation % description = 'mean percolation'
       sitesPercolation % unit = 'mm'
       sitesPercolation % offsetZ = 0.
       CALL WriteMetadata ( network = sitesPercolation, &
                        fileunit = fileUnitPercolation )
       CALL WriteData (sitesPercolation, fileUnitPercolation, .TRUE.)
   END IF
ELSE
   varOut (16) = .FALSE.
END IF

!actual evapotranspiration
IF ( IniReadInt ('actual-ET', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (et % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', 'et not allocated, &
                                        forced to not export basin average ')
       varOut (17) = .FALSE.
   ELSE
       varOut (17) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesET )
       fileUnitET = GetUnit ()
       string = TRIM(pathout) // 'mean_et.fts'
       OPEN ( unit = fileUnitET, file = string )
       sitesET % description = 'mean actual evapotranspiration'
       sitesET % unit = 'mm'
       sitesET % offsetZ = 0.
       CALL WriteMetadata ( network = sitesET, &
                        fileunit = fileUnitET )
       CALL WriteData (sitesET, fileUnitET, .TRUE.)
   END IF
ELSE
   varOut (17) = .FALSE.
END IF


!potential evapotranspiration
IF ( IniReadInt ('potential-ET', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (pet % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', 'pet not allocated, &
                                        forced to not export basin average ')
       varOut (18) = .FALSE.
   ELSE
       varOut (18) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesPET )
       fileUnitPET = GetUnit ()
       string = TRIM(pathout) // 'mean_pet.fts'
       OPEN ( unit = fileUnitPET, file = string )
       sitesPET % description = 'mean potential evapotranspiration'
       sitesPET % unit = 'mm'
       sitesPET % offsetZ = 0.
       CALL WriteMetadata ( network = sitesPET, &
                        fileunit = fileUnitPET )
       CALL WriteData (sitesPET, fileUnitPET, .TRUE.)
   END IF
ELSE
   varOut (18) = .FALSE.
END IF

!soil moisture variation (soil water storage variation)
IF ( IniReadInt ('delta-soil-moisture', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (dsm % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', 'delta-soil-moisture not allocated, &
                                        forced to not export basin average ')
       varOut (19) = .FALSE.
   ELSE
       varOut (19) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesDSM )
       fileUnitDSM = GetUnit ()
       string = TRIM(pathout) // 'mean_delta-soil-moisture.fts'
       OPEN ( unit = fileUnitDSM, file = string )
       sitesDSM % description = 'mean soil water storage variation'
       sitesDSM % unit = 'mm'
       sitesDSM % offsetZ = 0.
       CALL WriteMetadata ( network = sitesDSM, &
                        fileunit = fileUnitDSM )
       CALL WriteData (sitesDSM, fileUnitDSM, .TRUE.)
   END IF
ELSE
   varOut (19) = .FALSE.
END IF

!snow melt
IF ( IniReadInt ('snow-melt', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (snowmelt % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', 'snow-melt not allocated, &
                                        forced to not export basin average ')
       varOut (20) = .FALSE.
   ELSE
       varOut (20) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesSnowMelt )
       fileUnitSnowMelt = GetUnit ()
       string = TRIM(pathout) // 'mean_snow-melt.fts'
       OPEN ( unit = fileUnitSnowMelt, file = string )
       sitesSnowMelt % description = 'mean snow melt'
       sitesSnowMelt % unit = 'mm'
       sitesSnowMelt % offsetZ = 0.
       CALL WriteMetadata ( network = sitesSnowMelt, &
                        fileunit = fileUnitSnowMelt )
       CALL WriteData (sitesSnowMelt, fileUnitSnowMelt, .TRUE.)
   END IF
ELSE
   varOut (20) = .FALSE.
END IF


!ice water equivalent
IF ( IniReadInt ('ice-water-equivalent', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (iwe % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', 'ice-water-equivalent not allocated, &
                                        forced to not export basin average ')
       varOut (21) = .FALSE.
   ELSE
       varOut (21) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesIWE )
       fileUnitIWE = GetUnit ()
       string = TRIM(pathout) // 'mean_icewe.fts'
       OPEN ( unit = fileUnitIWE, file = string )
       sitesIWE % description = 'mean ice water equivalent'
       sitesIWE % unit = 'mm'
       sitesIWE % offsetZ = 0.
       CALL WriteMetadata ( network = sitesIWE, &
                        fileunit = fileUnitIWE )
       CALL WriteData (sitesIWE, fileUnitIWE, .TRUE.)
   END IF
ELSE
   varOut (21) = .FALSE.
END IF


!ice melt
IF ( IniReadInt ('ice-melt', iniDB ) == 1) THEN
   IF ( .NOT. ALLOCATED (icemelt % mat) ) THEN
       CALL Catch ('warning', 'BasinAverage', 'ice-melt not allocated, &
                                        forced to not export basin average ')
       varOut (22) = .FALSE.
   ELSE
       varOut (22) = .TRUE.
       countVar = countVar + 1
       CALL CopyNetwork ( sites, sitesIceMelt )
       fileUnitIceMelt = GetUnit ()
       string = TRIM(pathout) // 'mean_ice-melt.fts'
       OPEN ( unit = fileUnitIceMelt, file = string )
       sitesIceMelt % description = 'mean icemelt'
       sitesIceMelt % unit = 'mm'
       sitesIceMelt % offsetZ = 0.
       CALL WriteMetadata ( network = sitesIceMelt, &
                        fileunit = fileUnitIceMelt )
       CALL WriteData (sitesIceMelt, fileUnitIceMelt, .TRUE.)
   END IF
ELSE
   varOut (22) = .FALSE.
END IF





CALL IniClose (iniDB) 

!create grid to cumulate variable
IF ( countVar > 0 ) THEN
    CALL NewGrid (cum, mask)
END IF


RETURN
END SUBROUTINE InitBasinAverage




!==============================================================================
!| Description: 
!  Export basin averaged values
SUBROUTINE ExportBasinAverage   & 
!
 (time, temp, tmean, tmax, tmin, precipitation, &
  rh, radiation, netradiation, windspeed, daily_precipitation, &
  irrigation, swe, sm, runoff, infiltration, percolation, et, &
  pet, dsm, snowmelt, iwe, icemelt )  

IMPLICIT NONE

!arguments with intent in:  
TYPE (DateTime) , INTENT(IN) :: time
TYPE (grid_real), INTENT(IN) :: temp !!air temperarure (°C)
TYPE (grid_real), INTENT(IN) :: tmean !!air temperarure daily mean(°C)
TYPE (grid_real), INTENT(IN) :: tmax !!air temperarure daily max (°C)
TYPE (grid_real), INTENT(IN) :: tmin !!air temperarure daily min (°C)
TYPE (grid_real), INTENT(IN) :: precipitation !!precipitation rate (m/s)
TYPE (grid_real), INTENT(IN) :: rh !!air relative humidity (0-100)
TYPE (grid_real), INTENT(IN) :: radiation !!solar radiation (w/m2)
TYPE (grid_real), INTENT(IN) :: netradiation !!net radiation (w/m2)
TYPE (grid_real), INTENT(IN) :: windspeed !!wind speed (m/s)
TYPE (grid_real), INTENT(IN) :: daily_precipitation !!daily precipitation rate (m/s)
TYPE (grid_real), INTENT(IN) :: irrigation !!irrigation rate (m/s)
TYPE (grid_real), INTENT(IN) :: swe !!snow water equivalent (m)
TYPE (grid_real), INTENT(IN) :: sm !!soil mositure (-)
TYPE (grid_real), INTENT(IN) :: runoff !!runoff (m/s)
TYPE (grid_real), INTENT(IN) :: infiltration !!infiltration (m/s)
TYPE (grid_real), INTENT(IN) :: percolation !!deep percolation (m/s)
TYPE (grid_real), INTENT(IN) :: et !!actual evapotranspiration (m/s)
TYPE (grid_real), INTENT(IN) :: pet !!potential evapotranspiration (m/s)
TYPE (grid_real), INTENT(IN) :: dsm !!soil moisture variation (m)
TYPE (grid_real), INTENT(IN) :: snowmelt !!snow melt (m)
TYPE (grid_real), INTENT(IN) :: iwe !! ice water equivalent (m)
TYPE (grid_real), INTENT(IN) :: icemelt !! ice melt (m)
 

!local declarations
INTEGER (KIND = short) :: i, j

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

!precipitation
IF ( varOut (1) ) THEN
    !cumulate variable
    CALL BasinCumulate ( precipitation, 1000. * dtBasinAverage )
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesPrecipitation) 
    !set current time
    sitesPrecipitation % time = time
    !write data
    CALL WriteData (sitesPrecipitation, fileUnitPrecipitation ) 
END IF

!daily-precipitation
IF ( varOut (2) ) THEN
    !cumulate variable
    CALL BasinCumulate ( daily_precipitation, 1000. * 86400. )
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesPrecipitationDaily) 
    !set current time
    sitesPrecipitationDaily % time = time
    !write data
    CALL WriteData (sitesPrecipitationDaily, fileUnitPrecipitationDaily ) 
END IF

!air temperature
IF ( varOut (3) ) THEN
    !cumulate variable
    CALL BasinCumulate ( temp )
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesTemp) 
    !set current time
    sitesTemp % time = time
    !write data
    CALL WriteData (sitesTemp, fileUnitTemp ) 
END IF


!temperature-daily-mean
IF ( varOut (4) ) THEN
    !cumulate variable
    CALL BasinCumulate ( tmean )
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesTmean) 
    !set current time
    sitesTmean % time = time
    !write data
    CALL WriteData (sitesTmean, fileUnitTmean ) 
END IF

!temperature-daily-max
IF ( varOut (5) ) THEN
    !cumulate variable
    CALL BasinCumulate ( tmax )
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesTmax) 
    !set current time
    sitesTmax % time = time
    !write data
    CALL WriteData (sitesTmax, fileUnitTmax ) 
END IF


!temperature-daily-min
IF ( varOut (6) ) THEN
    !cumulate variable
    CALL BasinCumulate ( tmin )
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesTmin) 
    !set current time
    sitesTmin % time = time
    !write data
    CALL WriteData (sitesTmin, fileUnitTmin ) 
END IF

!relative-humidity
IF ( varOut (7) ) THEN
    !cumulate variable
    CALL BasinCumulate ( rh )
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesRH) 
    !set current time
    sitesRH % time = time
    !write data
    CALL WriteData (sitesRH, fileUnitRH ) 
END IF


!solar-radiation
IF ( varOut (8) ) THEN
    !cumulate variable
    CALL BasinCumulate ( radiation )
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesRadiation) 
    !set current time
    sitesRadiation % time = time
    !write data
    CALL WriteData (sitesRadiation, fileUnitRadiation ) 
END IF

!net-radiation
IF ( varOut (9) ) THEN
    !cumulate variable
    CALL BasinCumulate ( netradiation )
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesNetRadiation) 
    !set current time
    sitesNetRadiation % time = time
    !write data
    CALL WriteData (sitesNetRadiation, fileUnitNetRadiation ) 
END IF

!wind-speed
IF ( varOut (10) ) THEN
    !cumulate variable
    CALL BasinCumulate ( windspeed )
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesWindSpeed) 
    !set current time
    sitesWindSpeed % time = time
    !write data
    CALL WriteData (sitesWindSpeed, fileUnitWindSpeed ) 
END IF

!irrigation
IF ( varOut (11) ) THEN
    !cumulate variable
    CALL BasinCumulate ( irrigation )
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesIrrigation) 
    !set current time
    sitesIrrigation % time = time
    !write data
    CALL WriteData (sitesIrrigation, fileUnitIrrigation ) 
END IF

!snow water equivalent
IF ( varOut (12) ) THEN
    !cumulate variable
    CALL BasinCumulate ( swe, 1000. )
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesSWE) 
    !set current time
    sitesSWE % time = time
    !write data
    CALL WriteData (sitesSWE, fileUnitSWE ) 
    
END IF

!soil moisture
IF ( varOut (13) ) THEN
    !cumulate variable
    CALL BasinCumulate ( sm )
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesSM) 
    !set current time
    sitesSM % time = time
    !write data
    CALL WriteData (sitesSM, fileUnitSM ) 
END IF

!runoff
IF ( varOut (14) ) THEN
    !cumulate variable
    CALL BasinCumulate ( runoff, 1000. * dtBasinAverage )
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesRunoff) 
    !set current time
    sitesRunoff % time = time
    !write data
    CALL WriteData (sitesRunoff, fileUnitRunoff ) 
END IF

!infiltration
IF ( varOut (15) ) THEN
    !cumulate variable
    CALL BasinCumulate ( infiltration, 1000. * dtBasinAverage )
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesInfiltration) 
    !set current time
    sitesInfiltration % time = time
    !write data
    CALL WriteData (sitesInfiltration, fileUnitInfiltration ) 
END IF

!percolation
IF ( varOut (16) ) THEN
    !cumulate variable
    CALL BasinCumulate ( percolation, 1000. * dtBasinAverage )
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesPercolation) 
    !set current time
    sitesPercolation % time = time
    !write data
    CALL WriteData (sitesPercolation, fileUnitPercolation ) 
END IF

!actual evapotranspiration
IF ( varOut (17) ) THEN
    !cumulate variable
    CALL BasinCumulate ( et, 1000. * dtBasinAverage )
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesET) 
    !set current time
    sitesET % time = time
    !write data
    CALL WriteData (sitesET, fileUnitET ) 
END IF

!potential evapotranspiration
IF ( varOut (18) ) THEN
    !cumulate variable
    CALL BasinCumulate ( pet, 1000. * dtBasinAverage )
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesPET) 
    !set current time
    sitesPET % time = time
    !write data
    CALL WriteData (sitesPET, fileUnitPET ) 
END IF


!soil water storage variation
IF ( varOut (19) ) THEN
    !cumulate variable
    CALL BasinCumulate ( dsm, 1000.)
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesDSM) 
    !set current time
    sitesDSM % time = time
    !write data
    CALL WriteData (sitesDSM, fileUnitDSM ) 
END IF

!snow melt
IF ( varOut (20) ) THEN
    !cumulate variable
    CALL BasinCumulate ( snowmelt, 1000.)
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesSnowMelt) 
    !set current time
    sitesSnowMelt % time = time
    !write data
    CALL WriteData (sitesSnowMelt, fileUnitSnowMelt ) 
END IF

!ice water eqivalent
IF ( varOut (21) ) THEN
    !cumulate variable
    CALL BasinCumulate ( iwe, 1000.)
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesIWE) 
    !set current time
    sitesIWE % time = time
    !write data
    CALL WriteData (sitesIWE, fileUnitIWE ) 
END IF

!ice melt
IF ( varOut (22) ) THEN
    !cumulate variable
    CALL BasinCumulate ( icemelt, 1000.)
    !assign data to site network
    CALL AssignDataFromGrid ( cum, sitesIceMelt) 
    !set current time
    sitesIceMelt % time = time
    !write data
    CALL WriteData (sitesIceMelt, fileUnitIceMelt ) 
END IF

RETURN
END SUBROUTINE ExportBasinAverage
  
  
!==============================================================================
!| Description: 
!  cumulate variable along stream network
SUBROUTINE BasinCumulate   & 
!
 ( var, conv )  

IMPLICIT NONE

!arguments with intent in:    
TYPE (grid_real), INTENT (IN) :: var
REAL (KIND = float), OPTIONAL, INTENT (IN) :: conv !!conversion factor

!local declaration:
INTEGER (KIND = short) :: i, j, k, is, js
REAL (KIND = float) :: factor

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

!reset cumulated grid
cum = 0.

!cobversion factor
IF ( PRESENT (conv) ) THEN
    factor = conv
ELSE
    factor = 1.
END IF

DO k = 1, streamNetwork % nreach
    
  i = streamNetwork % branch (k) % i0
  j = streamNetwork % branch (k) % j0
  
  
  !follow the branch
  DO WHILE ( .NOT.( j == streamNetwork % branch (k) % j1 .AND. &
                    i == streamNetwork % branch (k) % i1  )    )
 
      
    !find downstream cell
    CALL DownstreamCell (i, j, flowDirection % mat(i,j), is, js )
   
    !if current cell is basin outlet exit
    IF ( CheckOutlet (i, j, is, js, flowDirection) ) THEN
        CYCLE
    END IF
    

    !add current cell
    cum % mat (i,j) = cum % mat (i,j) + var % mat (i,j) * factor
    
    !cumulate downstream
    cum % mat (is, js) = cum % mat (is, js) + cum % mat (i, j)
    
    !select downstream cell for next loop
    j = js
    i = is
  END DO  
END DO


!divide by number of cells in the basin to compute the mean
DO j = 1, cum % jdim
    DO i = 1, cum % idim
        IF ( cum % mat (i,j) /= cum % nodata ) THEN
            cum % mat (i,j) = cum % mat (i,j) / flowAccumulation % mat (i,j)
        END IF
    END DO
END DO

RETURN

 END SUBROUTINE BasinCumulate
 
 
 
 
!==============================================================================
!| Description: 
!  delineate and export basin mask
SUBROUTINE BasinRasterExport   & 
!
 ( pathout  )  

IMPLICIT NONE

!arguments with intent in:
CHARACTER (LEN = *), INTENT(IN) :: pathout


!local declarations
INTEGER (KIND = short) :: i
CHARACTER (LEN = 300) :: fileName
TYPE (grid_integer) :: basin
REAL (KIND = float) :: x, y

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

!initialize a dummy raster
CALL NewGrid (basin,flowDirection)

DO i = 1, sites % countObs
    !fileName = TRIM(pathout) // 'basin_mask_id_' // &
    fileName = 'basin_mask_id_' // TRIM (sites % obs (i) % id) // '.asc'
    x = sites % obs (i) % xyz % easting
    y = sites % obs (i) % xyz % northing
    CALL BasinDelineate (flowDirection,x,y, basin)
    CALL ExportGrid (basin, fileName, ESRI_ASCII)
END DO

CALL GridDestroy (basin)

RETURN
END SUBROUTINE BasinRasterExport

    
END MODULE BasinAverage