!! Manage a group of stations 
!|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.5 - 22nd July 2024  
!  
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 03/Jun/2011 | Original code |
! | 1.1      | 25/Nov/2016 | added WriteDataFileUnit routine |
! | 1.2      | 02/Feb/2023 | added AssignDataFromGrid |
! | 1.3      | 18/Feb/2023 | optional scaleFactor in AssignDataFromGrid |
! | 1.4      | 07/Jun/2024 | changed # with DateTime in output file header |
! | 1.5      | 22/Jul/2024 | updated the synchronization of time of data to read |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
!  set of fortran routines to manage a group of stations 
!  (surface meteorological or other stations) spread over 
!  a given area for making regular observations.  
! 
MODULE ObservationalNetworks


! Modules used: 
! 

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

USE Chronos, ONLY: &
  !imported definitions:
  DateTime, &
  !Imported operators:
  ASSIGNMENT( = ), &
  OPERATOR (<), &
  OPERATOR (>), &
  OPERATOR (>=), &
  OPERATOR (==), &
  !Imported variables:
  timeString, timeDefault

USE GridLib, ONLY : &
  !Imported definitions:
  grid_real

USE GridOperations, ONLY : &
  !Imported definitions:
  GetIJ

USE GeoLib, ONLY: &
  !Imported definition:
  Coordinate, CRS, &
  DecodeEPSG, &
  !Imported operators:
  ASSIGNMENT( = )

USE Utilities, ONLY: &
  !imported routines:
  GetUnit

USE StringManipulation, ONLY: &
  !Imported routines:
  StringToLower, ToString, &
  StringCompact

USE LogLib, ONLY: &
  !Imported routines:
  Catch

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

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

! Global Procedures:

PUBLIC :: ReadMetadata
PUBLIC :: WriteMetadata
PUBLIC :: ActualObservations
PUBLIC :: ReadData
PUBLIC :: WriteData
PUBLIC :: CopyNetwork
PUBLIC :: DestroyNetwork

! Global Type Definitions: 

! Global Declarations: 

! Local (i.e. private) Declarations: 
! Local Procedures:

PRIVATE :: ReadMetadataFileName
PRIVATE :: ReadMetadataFileUnit
PRIVATE :: WriteMetadataFileName
PRIVATE :: WriteMetadataFileUnit
PRIVATE :: ReadDataFileUnit
PRIVATE :: WriteDataFileUnit
PRIVATE :: AssignDataFromGridReal

! Operator definitions: 
!   Define new operators or overload existing ones.

INTERFACE ReadMetadata
   MODULE PROCEDURE  ReadMetadataFileName
   MODULE PROCEDURE  ReadMetadataFileUnit
END INTERFACE

INTERFACE WriteMetadata
   MODULE PROCEDURE  WriteMetadataFileName
   MODULE PROCEDURE  WriteMetadataFileUnit
END INTERFACE

INTERFACE ReadData
   MODULE PROCEDURE ReadDataFileUnit
END INTERFACE

INTERFACE WriteData
   MODULE PROCEDURE WriteDataFileUnit
END INTERFACE

INTERFACE AssignDataFromGrid
   MODULE PROCEDURE AssignDataFromGridReal
END INTERFACE

!!defintion of Observation type
TYPE Observation
    CHARACTER (LEN = 300)    :: name !!name of observational site
	REAL (KIND = float)      :: value !!The observation value itself
	CHARACTER (LEN = 20)     :: id !!identification code
    TYPE (Coordinate)        :: xyz !!easting, northing and elevation in real world
	!INTEGER (KIND = short)  :: col !!local easting coordinate in "raster world"
	!INTEGER (KIND = short)  :: row !!local northing coordinate in "raster world"
	REAL (KIND = float)      :: z   !!elevation [m] a.s.l.
END TYPE Observation

!!definition of ObservationalNetwork type
TYPE ObservationalNetwork
  TYPE(Observation), ALLOCATABLE :: obs (:) !!arbitrary observations in network
  INTEGER (KIND = short)     :: countObs !!number of observations in network
  REAL (KIND = float)        :: nodata !!conventional code for missing data
  INTEGER (KIND = short)     :: timeIncrement    ![s]
  TYPE (DateTime)            :: time !!The date and time of the observations
  INTEGER (KIND = short)     :: dataType !!1 = continuous, 2 = cumulative, 3 = incremental, 4 = average, 5 = maximum, 6 = minimum
  CHARACTER (LEN = 300)      :: description !!The name of the physical, chemical, or biological quantity that the value represents (e.g. streamflow, precipitation, water quality)
  CHARACTER (LEN = 100)      :: unit !!unit of measure
  LOGICAL                    :: syncData = .FALSE. !!true when network is synchronized to data section in file
  CHARACTER (LEN = 300)      :: path !!filesystem path of file containing observations
  REAL (KIND = float)        :: offsetZ !!offset from ground elevation [m]
  TYPE (CRS)                 :: mapping
  INTEGER (KIND = short)     :: epsg
END TYPE ObservationalNetwork

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

!==============================================================================
!| Description:
!   count number of observations different from missing data
!   Can optionally return network containing actual measurements.
SUBROUTINE ActualObservations &
!
(network, count, activeNetwork)

IMPLICIT NONE

!Arguments with intent (in):
TYPE (ObservationalNetwork), INTENT(IN) :: network

!Arguments with intent (out):
TYPE (ObservationalNetwork), OPTIONAL, INTENT(OUT) :: activeNetwork
INTEGER (KIND = short), INTENT (OUT) :: count

!local declarations
INTEGER (KIND = short) :: i
!--------------------------end of declarations---------------------------------
count = 0
DO i = 1, network % countObs
  IF (network % obs (i) % value /= network % nodata) THEN
    count = count + 1
  END IF
END DO

!active network
IF (PRESENT (activeNetwork)) THEN

    activeNetwork % countObs = count

    IF (activeNetwork % countObs == 0) THEN
      CALL Catch ('error', 'ObservationalNetworks', &
           'no available values to build active network' )   
    END IF

    activeNetwork % nodata = network % nodata 
    activeNetwork % timeIncrement = network % timeIncrement 
    activeNetwork % time = network % time 
    activeNetwork % dataType = network % dataType 
    activeNetwork % description = network % description 
    activeNetwork % unit = network % unit
    activeNetwork % offsetZ = network % offsetZ 
    activeNetwork % mapping = network % mapping

    ALLOCATE (activeNetwork % obs (count))

    !copy active values
    count = 0
    DO i = 1, network % countObs
      IF (network % obs (i) % value /= network % nodata) THEN
         count = count + 1
         activeNetwork % obs (count) % name = network % obs (i) % name
         activeNetwork % obs (count) % value = network % obs (i) % value
         activeNetwork % obs (count) % id = network % obs (i) % id
         activeNetwork % obs (count) % xyz = network % obs (i) % xyz
         activeNetwork % obs (count) % z = network % obs (i) % z
      END IF
    END DO		
END IF


END SUBROUTINE ActualObservations

!==============================================================================
!| Description:
!   make a copy of one observational network
SUBROUTINE CopyNetwork &
!
(input, output)

IMPLICIT NONE

! Arguments with intent (in):
TYPE (ObservationalNetwork), INTENT(IN) :: input

! Arguments with intent(inout):
TYPE (ObservationalNetwork), INTENT(INOUT) :: output

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

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

output % countObs      = input % countObs 
output % nodata        = input % nodata
output % timeIncrement = input % timeIncrement
output % time          = input % time
output % dataType      = input % dataType
output % description   = input % description
output % unit          = input % unit
output % mapping       = input % mapping
output % epsg          = input % epsg
output % offsetZ       = input % offsetZ

!allocate output network observations
IF (ALLOCATED (output % obs) ) THEN
   !check space is enough
   IF (output % countObs /= input % countObs) THEN
      DEALLOCATE (output % obs)
      ALLOCATE(output % obs (output % countObs))
   END IF
ELSE
    ALLOCATE(output % obs (output % countObs))
END IF

!create copy of observations
DO k = 1, output % countObs
  output % obs(k) % name         = input % obs(k) % name
  output % obs(k) % value        = input % obs(k) % value
  output % obs(k) % id           = input % obs(k) % id
  output % obs(k) % xyz          = input % obs(k) % xyz
!  output % obs(k) % col          = input % obs(k) % col
!  output % obs(k) % row          = input % obs(k) % row
  output % obs(k) % z            = input % obs(k) % z
END DO

RETURN
END SUBROUTINE CopyNetwork



!==============================================================================
!| Description:
!   read metadata section from file
SUBROUTINE ReadMetadataFileName &
!
(filename, network)

IMPLICIT NONE

! Arguments with intent (in):
CHARACTER (LEN = *), INTENT(IN) :: filename

! Arguments with intent(out):
TYPE (ObservationalNetwork), INTENT(OUT) :: network

!local declarations:
INTEGER (KIND = short) :: fileunit
INTEGER (KIND = short) :: err_io
INTEGER (KIND = short) :: i
CHARACTER (LEN = 300)  :: string
TYPE (IniList)         :: ini

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

!store path
network % path = filename


!open and read name value pairs name = value
CALL IniOpen (filename, ini)

!check and store network info
!description
IF (KeyIsPresent ('description', ini) ) THEN
    network % description = IniReadString ('description', ini)
ELSE !description is missing
    CALL Catch ('error', 'ObservationalNetworks',   &
                'description missing in file: ', &
                argument = TRIM(filename) )
END IF

!unit
IF (KeyIsPresent ('unit', ini) ) THEN
    network % unit = IniReadString ('unit', ini)
ELSE !unit is missing
    CALL Catch ('error', 'ObservationalNetworks',   &
                'unit missing in file: ', &
                argument = TRIM(filename) )
END IF

!epsg
IF (KeyIsPresent ('epsg', ini) ) THEN
    network % epsg = IniReadInt ('epsg', ini)
ELSE !epsg is missing
    CALL Catch ('error', 'ObservationalNetworks',   &
                'epsg missing in file: ', &
                argument = TRIM(filename) )
END IF

!assign coordinate reference system from EPSG code
network % mapping = DecodeEPSG (network % epsg)

!number of stations
IF (KeyIsPresent ('count', ini) ) THEN
    network % countObs = IniReadInt ('count', ini)
ELSE !count is missing
    CALL Catch ('error', 'ObservationalNetworks',   &
                'count missing in file: ', &
                argument = TRIM(filename) )
END IF

!dt
IF (KeyIsPresent ('dt', ini) ) THEN
    network % timeIncrement = IniReadInt ('dt', ini)
ELSE !dt is missing
    CALL Catch ('error', 'ObservationalNetworks',   &
                'dt missing in file: ', &
                argument = TRIM(filename) )
END IF

!missing-data flag
IF (KeyIsPresent ('missing-data', ini) ) THEN
    network % nodata = IniReadReal ('missing-data', ini)
ELSE !missing-data is missing
    CALL Catch ('error', 'ObservationalNetworks',   &
                'missing-data missing in file: ', &
                argument = TRIM(filename) )
END IF

!height offset
IF (KeyIsPresent ('offsetz', ini) ) THEN
    network % offsetZ = IniReadReal ('offsetz', ini)
ELSE !offsetz is missing
    CALL Catch ('error', 'ObservationalNetworks',   &
                'offsetz missing in file: ', &
                argument = TRIM(filename) )
END IF

!destroy ini
CALL IniClose (ini)

!open file
fileunit = GetUnit ()
OPEN (UNIT = fileunit, file = filename)

!scan for searching keyword "metadata"
string = ''
DO WHILE ( .NOT. (string(1:8) == 'metadata'))
!DO WHILE ( .NOT. (string(1:10) == 'anagrafica'))
	READ(fileunit,'(a)', iostat = err_io) string
	IF (err_io /= 0) THEN
	   CALL Catch ('error', 'ObservationalNetworks', &
       'keyword metadata not found in file: ',  &
       argument = TRIM(fileName) )     
    END IF
    string = StringCompact (string)
    string = StringToLower (string)
END DO

ALLOCATE (network % obs (network % countObs))

DO i = 1, network % countObs
	READ(fileunit,*) network % obs (i) % name, &
	                 network % obs (i) % id, &
	                 network % obs (i) % xyz % easting, &
	                 network % obs (i) % xyz % northing, &
	                 network % obs (i) % xyz % elevation
	!initialize observation values to nodata
	 network % obs (i) % value = network % nodata
END DO

CLOSE (fileunit)

RETURN
END SUBROUTINE ReadMetadataFileName


!==============================================================================
!| Description:
!   read metadata section from open unit
SUBROUTINE ReadMetadataFileUnit &
!
(fileunit, network)

IMPLICIT NONE

! Arguments with intent (in):
INTEGER (KIND = short), INTENT(IN) :: fileunit

! Arguments with intent(out):
TYPE (ObservationalNetwork), INTENT(OUT) :: network

!local declarations:
INTEGER (KIND = short) :: err_io
INTEGER (KIND = short) :: i
CHARACTER (LEN = 300)   :: string
TYPE (IniList)         :: ini
CHARACTER (LEN = 1000) :: filename
!----------------------------end of declarations-------------------------------

!store path
INQUIRE (UNIT = fileunit, NAME = network % path)

!open and read name value pairs name = value
CALL IniOpen (fileunit, ini)

!check and store network info
!description
IF (KeyIsPresent ('description', ini) ) THEN
    network % description = IniReadString ('description', ini)
ELSE !description is missing
    INQUIRE (fileunit, NAME = filename)
    CALL Catch ('error', 'ObservationalNetworks',   &
                'description missing in file: ', &
                argument = TRIM(filename) )
END IF

!unit
IF (KeyIsPresent ('unit', ini) ) THEN
    network % unit = IniReadString ('unit', ini)
ELSE !unit is missing
    INQUIRE (fileunit, NAME = filename)
    CALL Catch ('error', 'ObservationalNetworks',   &
                'unit missing in file: ', &
                argument = TRIM(filename) )
END IF

!epsg
IF (KeyIsPresent ('epsg', ini) ) THEN
    network % epsg = IniReadInt ('epsg', ini)
ELSE !epsg is missing
    INQUIRE (fileunit, NAME = filename)
    CALL Catch ('error', 'ObservationalNetworks',   &
                'epsg missing in file: ', &
                argument = TRIM(filename) )
END IF

!assign coordinate reference system from EPSG code
network % mapping = DecodeEPSG (network % epsg)

!number of stations
IF (KeyIsPresent ('count', ini) ) THEN
    network % countObs = IniReadInt ('count', ini)
ELSE !count is missing
    INQUIRE (fileunit, NAME = filename)
    CALL Catch ('error', 'ObservationalNetworks',   &
                'count missing in file: ', &
                argument = TRIM(filename) )
END IF

!dt
IF (KeyIsPresent ('dt', ini) ) THEN
    network % timeIncrement = IniReadInt ('dt', ini)
ELSE !dt is missing
    INQUIRE (fileunit, NAME = filename)
    CALL Catch ('error', 'ObservationalNetworks',   &
                'dt missing in file: ', &
                argument = TRIM(filename) )
END IF

!missing-data flag
IF (KeyIsPresent ('missing-data', ini) ) THEN
    network % nodata = IniReadReal ('missing-data', ini)
ELSE !missing-data is missing
    INQUIRE (fileunit, NAME = filename)
    CALL Catch ('error', 'ObservationalNetworks',   &
                'missing-data missing in file: ', &
                argument = TRIM(filename) )
END IF

!height offset
IF (KeyIsPresent ('offsetz', ini) ) THEN
    network % offsetZ = IniReadReal ('offsetz', ini)
ELSE !offsetz is missing
    INQUIRE (fileunit, NAME = filename)
    CALL Catch ('error', 'ObservationalNetworks',   &
                'offsetz missing in file: ', &
                argument = TRIM(filename) )
END IF

!destroy ini
CALL IniClose (ini)

!return to the beginning of file
REWIND (fileunit)

!scan for searching keyword "metadata"
string = ''
DO WHILE ( .NOT. (string(1:8) == 'metadata'))
!DO WHILE ( .NOT. (string(1:10) == 'anagrafica'))
	READ(fileunit,'(a)', iostat = err_io) string
	IF (err_io /= 0) THEN
	   CALL Catch ('error', 'ObservationalNetworks', &
       'keyword metadata not found' )     
    END IF
    string = StringCompact (string)
    string = StringToLower (string) 
END DO

ALLOCATE (network % obs (network % countObs))

DO i = 1, network % countObs
	READ(fileunit,*) network % obs (i) % name, &
	                 network % obs (i) % id, &
	                 network % obs (i) % xyz % easting, &
	                 network % obs (i) % xyz % northing, &
	                 network % obs (i) % xyz % elevation
	!initialize observation values to nodata
	 network % obs (i) % value = network % nodata
END DO

RETURN
END SUBROUTINE ReadMetadataFileUnit



!==============================================================================
!| Description:
!   write metadata section to file 
! TODO: full CRS support 
SUBROUTINE WriteMetadataFileName &
!
(network, filename)

IMPLICIT NONE

! Arguments with intent(in):
TYPE (ObservationalNetwork), INTENT(IN) :: network
CHARACTER (LEN = *), INTENT(IN) :: filename

!local declarations:
INTEGER (KIND = short) :: fileunit
INTEGER (KIND = short) :: i
!----------------------------end of declarations-------------------------------

fileunit = GetUnit ()

OPEN (UNIT = fileunit, file = filename)

WRITE(fileunit,'(a)') 'description = ', TRIM (network % description)
WRITE(fileunit,'(a)') 'unit = ', TRIM (network % unit)
WRITE(fileunit,'(a)') 'epsg = ', ToString (network % epsg)
WRITE(fileunit,'(a, i8)') 'count = ', network % countObs
WRITE(fileunit,'(a, i10)') 'dt = ', network % timeIncrement
WRITE(fileunit,'(a, f10.3)') 'missing-data = ', network % nodata
WRITE(fileunit,'(a, f10.3)') 'offsetz = ', network % offsetZ

WRITE(fileunit,*) !leave one blank line 
WRITE(fileunit,'(a)') 'metadata'
DO i = 1, network % countObs
	WRITE(fileunit,'(a30," ",a9," ",f12.2," ",f12.2," ",f12.2)') &
	      network % obs (i) % name, &
	      network % obs (i) % id,&
	      network % obs (i) % xyz % easting ,&
	      network % obs (i) % xyz % northing ,&
	      network % obs (i) % xyz % elevation
END DO

CLOSE (fileunit)

RETURN
END SUBROUTINE WriteMetadataFileName


!==============================================================================
!| Description:
!   write metadata section to open file unit 
! TODO: full CRS support 
SUBROUTINE WriteMetadataFileUnit &
!
(network, fileunit)

IMPLICIT NONE

! Arguments with intent(in):
TYPE (ObservationalNetwork), INTENT(IN) :: network
INTEGER (KIND = short), INTENT(IN) :: fileunit

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

WRITE(fileunit,'(a)') 'description = '// TRIM (network % description) 
WRITE(fileunit,'(a)') 'unit = ' // TRIM (network % unit)
WRITE(fileunit,'(a)') 'epsg = ' // ToString (network % epsg)
WRITE(fileunit,'(a, i8)') 'count = ', network % countObs
WRITE(fileunit,'(a, i10)') 'dt = ', network % timeIncrement
WRITE(fileunit,'(a, f10.3)') 'missing-data = ', network % nodata
WRITE(fileunit,'(a, f10.3)') 'offsetz = ', network % offsetZ

WRITE(fileunit,*) !leave one blank line 
WRITE(fileunit,'(a)') 'metadata'
DO i = 1, network % countObs
	WRITE(fileunit,'(a," ",a," ",f12.2," ",f12.2," ",f12.2)') &
	      TRIM ( network % obs (i) % name ), &
	      TRIM ( network % obs (i) % id ), &
	      network % obs (i) % xyz % easting ,&
	      network % obs (i) % xyz % northing ,&
	      network % obs (i) % xyz % elevation
END DO

RETURN
END SUBROUTINE WriteMetadataFileUnit



!==============================================================================
!| Description:
!   read data from file unit. Data spanned on multiple time steps
!   can be aggregated computing average, cumulated, maximum or minimum.
!   Aggregated value is considered as missing if number of 
!   actual available observations is less than a given percentage (tresh) 
SUBROUTINE ReadDataFileUnit &
!
(network, fileunit, time, aggr_time, aggr_type, tresh)

IMPLICIT NONE

! Arguments with intent (in):
INTEGER (KIND = short), INTENT(IN) :: fileunit

! Arguments with intent(inout):
TYPE (ObservationalNetwork), INTENT(INOUT) :: network

!Optional arguments:
TYPE (DateTime), OPTIONAL, INTENT(IN) :: time
INTEGER (KIND = short), OPTIONAL, INTENT(IN) :: aggr_time ![s]
CHARACTER (LEN = *), OPTIONAL, INTENT(IN) :: aggr_type !sum, ave, min, max
REAL (KIND = float), OPTIONAL, INTENT(IN) :: tresh ![0-1]

!local declarations:
CHARACTER (LEN = 10)   :: string
INTEGER (KIND = short) :: err_io
INTEGER (KIND = short) :: nstep
INTEGER (KIND = short) :: i, j
INTEGER (KIND = short) :: step
REAL (KIND = float)    :: valid
TYPE (ObservationalNetwork) :: tempNetwork
INTEGER (KIND = short), POINTER :: validObs(:)
CHARACTER (LEN = 300) :: filename
CHARACTER (LEN = 5) :: type

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

!Check optional arguments
IF( PRESENT(aggr_time) .AND. .not.PRESENT(aggr_type) ) THEN
    CALL Catch ('error', 'ObservationalNetworks', &
       'undefined aggregation type in routine ReadData')
END IF

!sync file to first data to read
!firstly scan for searching keyword "data"
IF ( .NOT. network % syncData) THEN
    string = ''
    DO WHILE ( .NOT. (string(1:4) == 'data'))
	    READ(fileunit,'(a)', iostat = err_io) string
	    IF (err_io /= 0) THEN
	       CALL Catch ('error', 'ObservationalNetworks', &
           'keyword data not found' )     
        END IF
        string = StringToLower (string)
    END DO
	READ(fileunit,*) !skip line
    
    !search for first data
    IF (PRESENT (time) ) THEN
        
        !scan file to search for time
        network % time = timeDefault
        DO WHILE (.not. (network % time >= time) )
            READ (fileunit,*,iostat = err_io) timeString
            
            network % time = timeString
            
        
            IF (err_io < 0) THEN !end of file reached

                INQUIRE (fileunit, NAME=filename)
			    timeString = time
	            CALL Catch ('error', 'ObservationalNetworks', &
                'time ' // timeString // ' not present in file: ', &
                argument = TRIM (filename) ) 
            END IF
        END DO    
        network % syncData = .TRUE.
        
    ELSE !sync file to first available data and return
        network % syncData = .TRUE.
        RETURN
    END IF
    
    !back to the preceding record.
    BACKSPACE (fileunit)
    IF ( network % time > time ) THEN ! back again
        BACKSPACE (fileunit)
    END IF
    
ELSE
   CALL Catch ('info', 'ObservationalNetworks', &
               ' network already synchronized to data section in file: ' , &
               argument = TRIM(network % path) ) 
END IF



IF( PRESENT(aggr_time) .AND. PRESENT(aggr_type) )THEN
        type = aggr_type
        type = StringToLower (type)
		IF((type /= 'ave').AND.(type /= 'sum').AND.&
		  (type /= 'min').AND.(type /= 'max'))THEN
		   CALL Catch ('error', 'ObservationalNetworks', &
            'wrong aggregation type in routine ReadData')
        END IF
        
        !initialize nstep
		IF( MOD(aggr_time,network % timeIncrement) == 0) THEN
			nstep = INT( aggr_time / network % timeIncrement )
		ELSE
		    CALL Catch ('error', 'ObservationalNetworks', &
            'wrong aggregation time in routine ReadData')
		END IF
END IF


!set threshold for missing data definition. Default value = 1 (100%)
IF (PRESENT(tresh)) THEN
	valid = tresh
ELSE
	valid = 1.
END IF

!set all observations value to missing data
DO i = 1, network % countObs
	network % obs (i) % value = network % nodata
END DO

!read data
IF( PRESENT(aggr_time) ) THEN
	CALL CopyNetwork (network, tempNetwork)
	ALLOCATE ( validObs(network % countObs) )
	validObs = nstep 
	DO step = 1, nstep
		READ(fileunit,*,iostat=err_io) timeString, &
		                (tempNetwork % obs(i) % value, i=1,tempNetwork % countObs)
		IF (err_io /= 0) THEN
	       CALL Catch ('error', 'ObservationalNetworks', &
           'error while reading data' )     
        END IF                
		network % time = timeString
		
		DO i = 1, tempNetwork % countObs
			IF (tempNetwork % obs (i) % value /= tempNetwork % nodata) THEN
				SELECT CASE (type(1:LEN_TRIM(type)))
					CASE ('ave','sum')
						IF (network % obs(i) % value == network % nodata) THEN
						    network % obs(i) % value = 0.
						END IF
						network % obs(i) % value = network % obs(i) % value + &
						                           tempNetwork % obs(i) % value
					CASE ('min')
						IF ((network % obs(i) % value == network % nodata) .OR. &
						   (network%obs(i)%value>tempNetwork%obs(i)%value))&
						    network%obs(i)%value=tempNetwork%obs(i)%value
					CASE ('max')
						IF ((network % obs(i) % value == network % nodata) .OR. &
						   (network % obs(i) % value < tempNetwork % obs(i) % value))&
						    network % obs(i) % value = tempNetwork % obs(i) % value
				END SELECT
			ELSE
				validObs(i) = validObs(i) - 1
			END IF
		END DO
	END DO
	
	SELECT CASE (type)
		CASE ('ave')
			DO i = 1, network % countObs
				IF((network % obs (i)% value /= network % nodata) .AND. &
				    (validObs(i) >= INT(valid*nstep)))     THEN
					network % obs(i) % value = network % obs(i) % value / validObs(i)
				ELSE
					network % obs(i) % value = network % nodata
				END IF
			END DO
		CASE ('sum')
			DO i = 1, network % countObs
				IF(.NOT.((network % obs(i) % value /= network % nodata) .AND. &
				    (validObs(i) >= INT(valid*nstep))))         THEN
					network % obs(i) % value = network % nodata
				END IF
			END DO
		CASE ('min','max')
			DO i = 1, network % countObs
				IF(.NOT.((network % obs(i) % value /= network % nodata) .AND. &
				         (validObs(i) >= INT (valid*nstep)))) THEN
					network % obs(i) % value = network % nodata
				END IF
			END DO
	END SELECT
	DEALLOCATE (tempNetwork % obs)
	DEALLOCATE (validObs)
ELSE
	READ(fileunit,*,iostat=err_io) timeString,(network % obs(j) % value,j = 1, network % countObs)
	IF (err_io /= 0) THEN
	   CALL Catch ('error', 'ObservationalNetworks', &
       'error reading  observation')     
    END IF
	
	network % time = timeString
	
END IF

RETURN
END SUBROUTINE ReadDataFileUnit




!==============================================================================
!| Description:
!   write data from file unit. if optional argument init is true
!   a string to initialize data section of file is written
SUBROUTINE WriteDataFileUnit &
!
(network, fileunit, init)

IMPLICIT NONE

! Arguments with intent (in):
TYPE (ObservationalNetwork), INTENT(IN) :: network
INTEGER (KIND = short), INTENT(IN) :: fileunit


!Optional arguments:
LOGICAL, OPTIONAL, INTENT(IN) :: init

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

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

!initialize data section
IF ( PRESENT (init) ) THEN
  IF (init) THEN
     WRITE (fileunit,*) !blank line
     WRITE (fileunit,'(a)') 'data' !data keyword
     WRITE (fileunit,'(a)',advance='no') 'DateTime'
				DO i = 1, network % countObs - 1
					WRITE (fileunit,'(" ",a)', advance='no') TRIM (network % obs (i) % id)
				END DO
        WRITE(fileunit,'(" ",a)') TRIM (network % obs (network % countObs) % id)
     RETURN !no need to go on
  END IF
END IF

!write current data
timeString = network % time
WRITE (fileunit,'(a25)',advance = 'no') timeString
DO i = 1, network % countObs - 1
	WRITE (fileunit,'(" ",e14.7)', advance = 'no') network % obs (i) % value
END DO
WRITE (fileunit,'(" ",e14.7)') network % obs (network % countObs) % value


RETURN
END SUBROUTINE WriteDataFileUnit



!==============================================================================
!| Description:
!   deallocate space 
SUBROUTINE DestroyNetwork &
!
(network)

IMPLICIT NONE

!Argument with intent (inout:
TYPE (ObservationalNetwork), INTENT (INOUT) :: network
!------------------------------------end of declarations-----------------------

DEALLOCATE (network % obs)

RETURN
END SUBROUTINE DestroyNetwork



!==============================================================================
!| Description:
!   assign data to a network from grid_real 
SUBROUTINE AssignDataFromGridReal &
!
( grid, network, scaleFactor  )

IMPLICIT NONE

!Arguments with intent(in)
TYPE (grid_real), INTENT (IN) :: grid
REAL (KIND = float), OPTIONAL, INTENT (IN) :: scaleFactor 

!Argument with intent (inout):
TYPE (ObservationalNetwork), INTENT (INOUT) :: network

!local declarations
INTEGER (KIND = short) :: i, j, k
REAL (KIND = float) :: x, y
LOGICAL             :: check
!------------------------------------end of declarations-----------------------

DO k = 1, network % countObs
    
    x = network % obs (k) % xyz % easting
    y = network % obs (k) % xyz % northing
    
    CALL GetIJ (x, y, grid, i, j, check)
    
    IF ( check .AND. grid % mat (i,j) /= grid % nodata ) THEN
        network % obs (k) % value = grid % mat (i,j)
    ELSE
        network % obs (k) % value =  network % nodata
    END IF
   
END DO

!apply scale factor
IF ( PRESENT ( scaleFactor ) ) THEN
    DO k = 1, network % countObs
      IF ( network % obs (k) % value /=  network % nodata ) THEN
           network % obs (k) % value = network % obs (k) % value * scaleFactor
      END IF
    END DO
END IF



RETURN
END SUBROUTINE AssignDataFromGridReal





END MODULE ObservationalNetworks