!! Simulate snow accumulation and melting 
!|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.3 - 15th January 2025 
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 15/Feb/2023 | Original code |
! | 1.1      | 29/May/2023 | change in configuration file to read initial snow water equivalent and water in snow |
! | 1.2      | 02/Dec/2024 | added refreezing of water retained in snow and snow water holding capacity |
! | 1.3      | 15/Jan/2025 | refreezing coefficient and snow water holding capacity read from configuration file |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! Module to simulate snow accumulation and melting.
! The code is a rewriting of the original code
! initially developed in the master thesis by
! Alessio Salandin and Giorgio Valè (2003)
!
MODULE Snow

! Modules used:

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

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

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

USE LogLib, ONLY: &
  !Imported routines:
  Catch

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

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

USE AirTemperature, ONLY: &
    !imported variables:
    temperatureAir

USE Precipitation, ONLY: &
   !imported variables:
   precipitationRate

USE StringManipulation, ONLY : &
    !Imported routines:
    ToString

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

USE MorphologicalProperties, ONLY: &
  !imported variables:
  flowDirection, dem

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

USE Utilities, ONLY : &
  !Imported routines:
  GetUnit

!Global variables:
INTEGER (KIND = short) :: dtSnow  !!computation time step
TYPE (grid_real) :: swe !! snow water equivalent (m)
TYPE (grid_real) :: rainfall !! liquid precipitation (m/s)
TYPE (grid_real) :: meltCoefficient !! melt coefficient (mm/day/°C) 
TYPE (grid_real) :: meltTemperature !! threshold temperature for melt starts (°C) 
TYPE (grid_real) :: upperTemperature !! upper temperature for partitioning precipitation (°C)
TYPE (grid_real) :: lowerTemperature !! lower temperature for partitioning precipitation (°C)
TYPE (grid_real) :: snowConductivity !! snow hydraulic conductivity (m/s)
TYPE (grid_real) :: waterInSnow  !! free water in snow pack (m)
TYPE (grid_real) :: swhc !!snow water holding capacity (-)
TYPE (grid_real) :: cRefreeze  !!coefficient of refreeze (-)
TYPE (grid_real) :: excessWaterSnowFlux  !! excess of water retained in snow flux (m/s)
TYPE (grid_real) :: snowMelt   !! snow melt amount within the time step (m)
TYPE (grid_real) :: rainfallrate   !!  precipitation liquid fraction (m/s)

!Global routines:
PUBLIC :: SnowInit
PUBLIC :: SnowUpdate
PUBLIC :: SnowPointInit
PUBLIC :: DegreeDay

!local variables:
TYPE (grid_integer)        , PRIVATE :: meltModel  !!melt model map
TYPE (grid_real)           , PRIVATE :: QinSnow, QoutSnow
TYPE (ObservationalNetwork), PRIVATE :: sites
TYPE (ObservationalNetwork), PRIVATE :: sitesSWE
INTEGER (KIND = short)     , PRIVATE :: fileUnitPointSWE
TYPE (DateTime)            , PRIVATE :: timePointExport


!local routines
PRIVATE :: SnowParameterUpdate
PRIVATE :: SnowPointExport
PRIVATE :: Refreezing


!=======
CONTAINS
!=======
! Define procedures contained in this module. 
  
!==============================================================================
!| Description: 
!  Initialize snow model
SUBROUTINE SnowInit   & 
  !
  (inifile, mask, time)       

IMPLICIT NONE

!Arguments with intent(in):
CHARACTER (LEN = *), INTENT(IN) :: inifile !!stores configuration information
TYPE (grid_integer), INTENT(IN) :: mask
TYPE (DateTime),     INTENT(IN) :: time

!local declarations:
TYPE (IniList)         :: iniDB
REAL (KIND = float)    :: scalar
INTEGER (KIND = short) :: iscalar

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

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

!load melt model
IF (SectionIsPresent('melt-model', iniDB)) THEN	
    IF (KeyIsPresent ('scalar', iniDB, 'melt-model') ) THEN
        iscalar = IniReadInt ('scalar', iniDB, 'melt-model')
        CALL NewGrid (meltModel, mask, iscalar)
    ELSE
         CALL GridByIni (iniDB, meltModel, section = 'melt-model')
    END IF
   
ELSE
    CALL Catch ('error', 'Snow',   &
			      'melt-model not found in configuration file' )
END IF

!set melt coefficient
IF (SectionIsPresent('melt-coefficient', iniDB)) THEN	
    IF (KeyIsPresent ('scalar', iniDB, 'melt-coefficient') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'melt-coefficient')
        CALL NewGrid ( meltCoefficient, mask, scalar )
    ELSE
         CALL GridByIni (iniDB, meltCoefficient, &
                         section = 'melt-coefficient', &
                         time = time )
    END IF
   
ELSE
    CALL Catch ('error', 'Snow',   &
			      'melt-coefficient not found in configuration file' )
END IF

!set threshold temperature for melt starts (°C)
IF (SectionIsPresent('melt-threshold-temperature', iniDB)) THEN	
    IF (KeyIsPresent ('scalar', iniDB, 'melt-threshold-temperature') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'melt-threshold-temperature')
        CALL NewGrid ( meltTemperature, mask, scalar )
    ELSE
         CALL GridByIni (iniDB, meltTemperature, &
                         section = 'melt-threshold-temperature', &
                         time = time )
    END IF
   
ELSE
    CALL Catch ('error', 'Snow',   &
        'melt-threshold-temperature not found in configuration file' )
END IF


!set upper temperature for partitioning precipitation (°C)
IF (SectionIsPresent('partitioning-upper-temperature', iniDB)) THEN	
    IF (KeyIsPresent ('scalar', iniDB, 'partitioning-upper-temperature') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'partitioning-upper-temperature')
        CALL NewGrid ( upperTemperature, mask, scalar )
    ELSE
        CALL GridByIni (iniDB, upperTemperature, &
                         section = 'partitioning-upper-temperature', &
                         time = time )
    END IF
   
ELSE
    CALL Catch ('error', 'Snow',   &
			      'partitioning-upper-temperature not found in configuration file' )
END IF


!set lower temperature for partitioning precipitation (°C)
IF (SectionIsPresent('partitioning-lower-temperature', iniDB)) THEN	
    IF (KeyIsPresent ('scalar', iniDB, 'partitioning-lower-temperature') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'partitioning-lower-temperature')
        CALL NewGrid( lowerTemperature, mask, scalar )
    ELSE
         CALL GridByIni (iniDB, lowerTemperature, &
                         section = 'partitioning-lower-temperature', &
                         time = time )
    END IF
   
ELSE
    CALL Catch ('error', 'Snow',   &
			      'partitioning-lower-temperature not found in configuration file' )
END IF

!set snow hydraulic conductivity (m/s)
IF (SectionIsPresent('hydraulic-conductivity', iniDB)) THEN	
    IF (KeyIsPresent ('scalar', iniDB, 'hydraulic-conductivity') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'hydraulic-conductivity')
        CALL NewGrid ( snowConductivity, mask, scalar )
    ELSE
         CALL GridByIni ( iniDB, snowConductivity, &
                         section = 'hydraulic-conductivity', &
                         time = time )
    END IF
   
ELSE
    CALL Catch ('error', 'Snow',   &
			      'hydraulic-conductivity not found in configuration file' )
END IF

!set refreezing coefficient (-)
IF (SectionIsPresent('refreezing-coefficient', iniDB)) THEN	
    IF (KeyIsPresent ('scalar', iniDB, 'refreezing-coefficient') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'refreezing-coefficient')
        CALL NewGrid ( cRefreeze, mask, scalar )
    ELSE
         CALL GridByIni ( iniDB, cRefreeze, &
                         section = 'refreezing-coefficient' )
    END IF
   
ELSE !set default value
    CALL NewGrid ( cRefreeze, mask, 0.05 )
END IF


!set snow water holding capacity (-)
IF (SectionIsPresent('water-holding-capacity', iniDB)) THEN	
    IF (KeyIsPresent ('scalar', iniDB, 'water-holding-capacity') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'water-holding-capacity')
        CALL NewGrid ( swhc, mask, scalar )
    ELSE
         CALL GridByIni ( iniDB, swhc, &
                         section = 'water-holding-capacity' )
    END IF
   
ELSE !set default value
    CALL NewGrid ( swhc, mask, 0.1 )
END IF


!set initial optional variables

!swe
IF (SectionIsPresent('snow-water-equivalent', iniDB)) THEN	
    IF (KeyIsPresent ('scalar', iniDB, 'snow-water-equivalent') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'snow-water-equivalent')
        CALL NewGrid ( swe, mask, scalar )
    ELSE
         CALL GridByIni (iniDB, swe, section = 'snow-water-equivalent')
    END IF
   
ELSE !set to default = 0
   CALL NewGrid ( swe, mask, 0. )
END IF

!water in snow
IF (SectionIsPresent('water-in-snow', iniDB)) THEN	
    IF (KeyIsPresent ('scalar', iniDB, 'water-in-snow') ) THEN
        scalar = IniReadReal ('scalar', iniDB, 'water-in-snow')
        CALL NewGrid ( waterInSnow, mask, scalar )
    ELSE
         CALL GridByIni (iniDB, waterInSnow, section = 'water-in-snow')
    END IF
   
ELSE !set to default = 0
   CALL NewGrid ( waterInSnow, mask, 0. )
END IF

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

!allocate variables
CALL NewGrid ( excessWaterSnowFlux, mask, 0. )
CALL NewGrid ( snowMelt,    mask, 0. )
CALL NewGrid ( QinSnow,     mask, 0. )
CALL NewGrid ( QoutSnow,    mask, 0. )




RETURN
END SUBROUTINE SnowInit
  
  
  
!==============================================================================
!| Description: 
!  compute accumulation and melting and update snow water equivalent
SUBROUTINE SnowUpdate   & 
  !
  ( time, mask )       

IMPLICIT NONE

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

!local declarations:
INTEGER (KIND = short) :: i, j, is, js
REAL (KIND = float) :: meltRate  !!snow melt rate (m/s)
REAL (KIND = float) :: refreezingRate !!refreezing rate of water retained in snowpack (m/s)
!REAL (KIND = float) :: cRefreeze = 0.05!!coefficient of refreeze (-)
!REAL (KIND = float) :: swhc = 0.1 !!snow water holding capacity (-)
REAL (KIND = float) :: alpha !!precipitation partitioning factor
REAL (KIND = float) :: liquid !!liquid fraction of precipitation (rain) (m/s)
REAL (KIND = float) :: solid !!solid fraction of precipitation (snow) (m/s)
REAL (KIND = float) :: tlower
REAL (KIND = float) :: tupper
REAL (KIND = float) :: tAir
REAL (KIND = float) :: cmelt !!snow melt coefficient (m/s/°C)
REAL (KIND = float) :: tmelt
REAL (KIND = float) :: swe_tminus1 !!swe at previous time step
REAL (KIND = float) :: ddx	!!actual cell length (m)
REAL (KIND = float) :: ds !!thickness of the saturated depth (m)
REAL (KIND = float) :: cellWidth !!cell width (m)
REAL (KIND = float) :: ijSlope  !!local slope (m/m)
REAL (KIND = float) :: qin, qout  !!input and output water discharge in snowpack
REAL (KIND = float) :: area !!cell area (m2)
!------------------------------end of declarations-----------------------------  

!check if parameter update is required
CALL SnowParameterUpdate (time)

!initialize melt grid
snowMelt = 0.


!computer inter-cell lateral water flux
QinSnow = 0.
QoutSnow = 0.
DO i = 1, mask % idim
    DO j = 1, mask % jdim
        IF ( mask % mat (i,j) == mask % nodata ) CYCLE
        
        !set downstream cell  (is,js)
        CALL DownstreamCell ( i, j, flowDirection % mat(i,j), is, js, ddx, &
                            flowDirection ) 
        
        IF ( CheckOutlet (i, j, is, js, flowDirection) ) CYCLE
      
        !thickness of the saturated depth
        ds = waterInSnow % mat (i,j)
        
        !cell width
        cellWidth = ( CellArea (mask, i, j) ) ** 0.5
        
        !local slope
        ijSlope = ( dem % mat (i,j) - dem % mat (is,js) ) / ddx
        
        IF ( ijSlope <= 0. ) THEN
            ijSlope = 0.0001
        END IF
        
        QoutSnow % mat (i,j) = cellWidth * ds * &
                               snowConductivity % mat (i,j) *  ijSlope
      
        !output becomes input of downstream cell
        QinSnow % mat (is,js) = QinSnow % mat (is,js) + QoutSnow % mat (i,j)
        
    END DO
END DO

!loop across cells
DO i = 1, mask % idim
    DO j = 1, mask % jdim
        
        !skip nodata
        IF ( mask % mat (i,j) == mask % nodata ) THEN
            CYCLE
        END IF
        
        !potential snow melt rate
        tAir  = temperatureAir % mat (i,j)
        cmelt = meltCoefficient % mat (i,j) / 1000. / 86400.
        tmelt = meltTemperature % mat (i,j)
        SELECT CASE ( meltModel % mat (i,j) )
        CASE (1) !degree-day temperature based
            meltRate = DegreeDay ( tAir,  tmelt,  cmelt )
        CASE DEFAULT
            CALL Catch ('error', 'Snow', 'melt model not implemented: ', &
                         argument = ToString (meltModel % mat (i,j) ) )
        END SELECT
        
        ! actual snow melt, limited by available snow 
        ! in the previous time step
        meltRate = MIN ( meltRate, swe % mat (i,j) / dtSnow )
        snowMelt % mat (i,j) = meltRate * dtSnow
        
        !precipitation partitioning 
        tupper = upperTemperature % mat (i,j)
        tlower = lowerTemperature % mat (i,j)
        
        IF ( tAir <= tlower ) THEN
            alpha = 0.
        ELSE IF ( tAir > tlower .AND.  tAir < tupper ) THEN
            alpha = ( tAir - tlower ) / ( tupper - tlower )
		ELSE
			alpha  = 1.
        END IF
        
        liquid = alpha * precipitationRate % mat (i,j)
        solid  = ( 1. - alpha ) * precipitationRate % mat (i,j)
        
        ! potential refreezing rate
        refreezingRate = Refreezing ( tAir,  tmelt,  cmelt, cRefreeze % mat (i,j) )
        
         ! actual refreezing, limited by available retained water   
        ! in the previous time step
        refreezingRate = MIN ( refreezingRate, &
                               waterInSnow % mat (i,j) / dtSnow )
        
        !update snow water equivalent
        swe % mat (i,j) = swe % mat (i,j) + &
                       ( solid - meltRate + refreezingRate ) *  dtSnow
        
        !update water in snow
        area = CellArea (mask, i, j) 
        qin = QinSnow % mat (i,j) / area
        qout = QoutSnow % mat (i,j) / area
        IF ( swe % mat (i,j) > 0. ) THEN
            
            waterInSnow % mat (i,j) = waterInSnow % mat (i,j) + &
                                      snowMelt % mat (i,j) +  &
                                      liquid * dtSnow + &
                                      ( qin - qout ) * dtSnow + &
                                     -  refreezingRate  * dtSnow 
            
            
            rainfallRate % mat (i,j) = 0.
            
        ELSE !rainfall does not contribute 
            
            waterInSnow % mat (i,j) = waterInSnow % mat (i,j) + &
                                      snowMelt % mat (i,j) +  &
                                      ( qin - qout ) * dtSnow  + &
                                      -  refreezingRate  * dtSnow
            
             rainfallRate % mat (i,j) = liquid
             
        END IF
        
        IF ( waterInSnow % mat (i,j) < 0. ) THEN
            waterInSnow % mat (i,j) = 0.
        END IF
        
        ! check retained water does not exceed water holding capacity of snow
        IF ( waterInSnow % mat (i,j) > swhc % mat (i,j) * swe % mat (i,j) ) THEN
            !exceeding water becomes runoff
            excessWaterSnowFlux % mat (i,j) = ( waterInSnow % mat (i,j) - &
                                            swhc % mat (i,j) * swe % mat (i,j) ) / dtSnow
            !set maximum 
            waterInSnow % mat (i,j) = swhc % mat (i,j) * swe % mat (i,j)
            
        END IF
        
        
    END DO
END DO

!write point SWE
IF ( time == timePointExport ) THEN
   CALL SnowPointExport ( time )
END IF

RETURN
END SUBROUTINE SnowUpdate
  
  
  
!==============================================================================
!| Description: 
!  update parameter map that change in time
SUBROUTINE SnowParameterUpdate   & 
  !
  (time)       

IMPLICIT NONE

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

!local declarations:
CHARACTER (LEN = 300) :: filename
CHARACTER (LEN = 300) :: varname

!------------------------------end of declarations-----------------------------
  
  
  
!melt coefficient
IF (  time == meltCoefficient % next_time ) THEN
   !destroy current grid
   filename = meltCoefficient % file_name
   varname = meltCoefficient % var_name
   CALL GridDestroy (meltCoefficient)
   !read new grid in netcdf file
   CALL NewGrid (meltCoefficient, TRIM(filename), NET_CDF, &
                  variable = TRIM(varname), time = time)
END IF

!melt temperature
IF (  time == meltTemperature % next_time ) THEN
   !destroy current grid
   filename = meltTemperature % file_name
   varname = meltTemperature % var_name
   CALL GridDestroy (meltTemperature)
   !read new grid in netcdf file
   CALL NewGrid (meltTemperature, TRIM(filename), NET_CDF, &
                  variable = TRIM(varname), time = time)
END IF

!upper temperature
IF (  time == upperTemperature % next_time ) THEN
   !destroy current grid
   filename = upperTemperature % file_name
   varname = upperTemperature % var_name
   CALL GridDestroy (upperTemperature)
   !read new grid in netcdf file
   CALL NewGrid (upperTemperature, TRIM(filename), NET_CDF, &
                  variable = TRIM(varname), time = time)
END IF

!lower temperature
IF (  time == lowerTemperature % next_time ) THEN
   !destroy current grid
   filename = lowerTemperature % file_name
   varname = lowerTemperature % var_name
   CALL GridDestroy (lowerTemperature)
   !read new grid in netcdf file
   CALL NewGrid (lowerTemperature, TRIM(filename), NET_CDF, &
                  variable = TRIM(varname), time = time)
END IF

!hydraulic conductivity
IF (  time == snowConductivity % next_time ) THEN
   !destroy current grid
   filename = snowConductivity % file_name
   varname = snowConductivity % var_name
   CALL GridDestroy (snowConductivity)
   !read new grid in netcdf file
   CALL NewGrid (snowConductivity, TRIM(filename), NET_CDF, &
                  variable = TRIM(varname), time = time)
END IF
  
RETURN
END SUBROUTINE SnowParameterUpdate

  
!==============================================================================
!| Description: 
!  compute melt rate (m/s) using degreeday temperature based model
FUNCTION DegreeDay   & 
  !
  ( tempAir, tempMelt, cMelt ) &
  !
RESULT (melt)

IMPLICIT NONE

!Arguments with intent(in):
REAL (KIND = float), INTENT(IN) :: tempAir !! air temperature (°C)
REAL (KIND = float), INTENT(IN) :: tempMelt !! melting temperature (°C)
REAL (KIND = float), INTENT(IN) :: cMelt !! melting coefficient (m/s/°C)

!local declarations:
REAL (KIND = float) :: melt !! melt rate (m/s)

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

IF ( tempAir <= tempMelt ) THEN
    melt = 0.
ELSE
   melt = cMelt * ( tempAir - tempMelt )	
END IF


RETURN
END FUNCTION DegreeDay


!==============================================================================
!| Description: 
!  compute refreezing rate (m/s) of the water retained in snowpack
!
! References:
!   van Verseveld, W. J., Weerts, A. H., Visser, M., Buitink, J., 
!   Imhoff, R. O., Boisgontier, H., Bouaziz, L., Eilander, D., Hegnauer, M., 
!   ten Velden, C., and Russell, B.: Wflow_sbm v0.7.3, a spatially distributed 
!   hydrological model: from global data to local applications, Geosci. 
!   Model Dev., 17, 3199–3234, https://doi.org/10.5194/gmd-17-3199-2024, 2024
!
FUNCTION Refreezing   & 
  !
  ( tempAir, tempMelt, cMelt, cRefreeze ) &
  !
RESULT (freeze)

IMPLICIT NONE

!Arguments with intent(in):
REAL (KIND = float), INTENT(IN) :: tempAir !! air temperature (°C)
REAL (KIND = float), INTENT(IN) :: tempMelt !! melting temperature (°C)
REAL (KIND = float), INTENT(IN) :: cMelt !! melting coefficient (m/s/°C)
REAL (KIND = float), INTENT(IN) :: cRefreeze !! refreezing coefficient (-)

!local declarations:
REAL (KIND = float) :: freeze !! refreezing rate (m/s)

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

IF ( tempAir <= tempMelt ) THEN
    freeze = cRefreeze * cMelt * ( tempMelt - tempAir )	   
ELSE
   freeze = 0.
END IF


RETURN
END FUNCTION Refreezing



!==============================================================================
!| Description:
!   Initialize export of point site data
SUBROUTINE SnowPointInit &
!
( pointfile, path_out, time )

IMPLICIT NONE

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

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

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

timePointExport = time

!open point file

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

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

!Read metadata
CALL ReadMetadata (fileunit, sites)

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

CLOSE ( fileunit )

!create virtual network and initialize file for output
fileUnitPointSWE = GetUnit ()
OPEN ( unit = fileUnitPointSWE, &
    file = TRIM(path_out) // 'point_swe.fts' )
    
CALL CopyNetwork ( sites, sitesSWE )

sitesSWE % description = 'snow water equivalent data exported from FEST'

sitesSWE % unit = 'mm'

sitesSWE % offsetZ = 0.

CALL WriteMetadata ( network = sitesSWE, &
                fileunit = fileUnitPointSWE )

CALL WriteData (sitesSWE, fileUnitPointSWE, .TRUE.)    

! destroy sites
CALL DestroyNetwork ( sites )
RETURN
END SUBROUTINE SnowPointInit

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

IMPLICIT NONE

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

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


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



!set current time
sitesSWE % time = time 
    
!populate data
CALL AssignDataFromGrid ( swe, sitesSWE, scaleFactor = 1000. )
    
!write data
CALL WriteData ( sitesSWE, fileUnitPointSWE )
    
timePointExport = timePointExport + sitesSWE % timeIncrement


RETURN
END SUBROUTINE SnowPointExport


END MODULE Snow