!! Manage mean daily air temperature
!|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.0 - 19th May 2021  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 19/May/2021 | Original code |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! Module to manage mean daily air temperature
!
MODULE AirTemperatureDailyMean
 

! Modules used: 

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


USE Loglib, ONLY : &
  !Imported routines:
  Catch

USE ObservationalNetworks, ONLY : &
  !Imported types:
  ObservationalNetwork, &
  !Imported routines:
  ReadMetadata, ReadData, &
  CopyNetwork, &
  !Imported operators:
  ASSIGNMENT (=)

USE Chronos, ONLY : &
  !Imported types:
  DateTime, &
  !Imported variables:
  timeString, &
  !Imported operations:
  OPERATOR (<), &
  OPERATOR (>), &
  OPERATOR (>=), &
  OPERATOR (<=), &
  OPERATOR (+), &
  OPERATOR (==), &
  ASSIGNMENT (=)

USE GridLib, ONLY : &
  !Imported types:
  grid_real, grid_integer, &
  !Imported routines:
  NewGrid, GetDtGrid, &
  GridDestroy, ExportGrid, &
  SetLongName, SetStandardName, &
  SetUnits, SetCurrentTime, &
  !Imported parameters:
  NET_CDF, ESRI_ASCII, &
  ESRI_BINARY
  

USE GridOperations, ONLY : &
  !Imported routines:
  GridByIni, GridResample, &
  CRSisEqual, GridConvert, &
  !Imported operations:
  ASSIGNMENT (=)

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

USE Utilities, ONLY: &
  !Imported routines:
  GetUnit

USE MeteoUtilities, ONLY: &
  !Imported routines:
  ReadField, ShiftMeteoWithLapse, &
  Interpolate, ShiftBackWithLapse, &
  Offset, Scale

USE FileSys, ONLY: &
  !Imported routines:
  DirExists, FileExists, &
  FileDelete, CurrentDir, &
  GetOS, DirNew, &
  !Imported parameters:
  WIN32

USE GeoLib, ONLY: &
  !Imported variables:
  point1, point2, &
  !Imported routines:
  DecodeEPSG, Convert, &
  !Imported types:
  CRS, &
  !Imported operations:
  OPERATOR (==)

USE StringManipulation, ONLY: &
   !Imported routines:
   ReplaceChar, ToString

USE Statistics, ONLY: &
    !imported routines:
    LinearRegressionSlope

USE Units, ONLY: &
    !imported parameters:
    day
    
IMPLICIT NONE 

! Global (i.e. public) Declarations: 

INTEGER (KIND = short)     :: dtTemperatureDailyMean  = 0 !!cumulation deltat 
TYPE(ObservationalNetwork) :: thermometersDailyMean !!temperature stations network
TYPE(grid_real)            :: temperatureAirDailyMean	   	  !![°C]


!Private declarations
TYPE(DateTime), PRIVATE         :: timeNew !!time when new  data must be read
TYPE(ObservationalNetwork), PRIVATE :: stationsRefElev !!stations network at reference elevation
INTEGER (KIND = short), PRIVATE :: fileunit !! unit of file containing  data
INTEGER (KIND = short), PRIVATE :: interpolationMethod  !!method to spatial interpolate site data
INTEGER (KIND = short), PRIVATE :: interpolationMethod_assignment  !!method to assign spatial interpolation !
                                                                   !!1 = one method for the entire domain, 2 = a map with interpolation method codes
TYPE (grid_integer),    PRIVATE :: interpolationMethod_map
INTEGER (KIND = short), PRIVATE :: interpolationMethod_vector (3) !!defines active interpolation methods 
TYPE (grid_real),       PRIVATE :: interpolatedMap (3)  !!1 = map for thiessen, 2 = map for idw, 3 = map for kriging
REAL (KIND = float),    PRIVATE :: cellsizeInterpolation !!spatial resolution of interpolated grid
TYPE(grid_real),        PRIVATE :: gridTemp !temporary grid used to interpolate to 
                                                     !different spatial resolution
INTEGER (KIND = short), PRIVATE :: neighbors  !!number of closest data to use for interpolation
TYPE(grid_real),        PRIVATE :: grid_devst !!standard deviation of kriging interpolation
CHARACTER (LEN = 300),  PRIVATE :: fileGrid !!file containing grids
INTEGER (KIND = short), PRIVATE :: dtGrid !!dt of imported  field
INTEGER (KIND = short), PRIVATE :: elevationDrift !! 1 = use elevation to modify interpolated data
TYPE(grid_real),        PRIVATE :: lapse_map !!lapse rate map
REAL (KIND = float),    PRIVATE :: lapse  !!lapse rate [°C/m]
REAL (KIND = float),    PRIVATE :: lapse_computed  !!lapse rate computed from data [°C/m]
REAL (KIND = float),    PRIVATE :: r2min  !! value of linear regression R2 below which lapse-rate-scalar is used instead of computed lapse rate
INTEGER (KIND = short), PRIVATE :: lapse_rate_assignment !!1 = scalar, 2 = map that may change with time
INTEGER (KIND = short), PRIVATE :: lapse_rate_computation !!1 = compute lapse rate from data at each time step, 0 = use assigned lapse rate
REAL, ALLOCATABLE, PRIVATE      :: vectorT (:), vectorZ (:) !!vector of temperatures and elevation to compute lapse rate
INTEGER (KIND = short), PRIVATE :: i, j  !!loop index
REAL (KIND = float), PARAMETER, PRIVATE :: refElevation = 1000. !!reference elevation for applying lapse rate
REAL (KIND = float), PRIVATE, PARAMETER :: MISSING_DEF_REAL = -9999.9
REAL (KIND = float), PRIVATE  :: scale_factor  !!scale factor to apply to final map
REAL (KIND = float), PRIVATE  :: offset_value !!offset to apply to final map
REAL (KIND = float), PRIVATE :: valid_prcn !!when data from several time steps are read, 
                                      !!this is the minimum percentage (0-1) of valid 
                                       !!data that must be prresent to consider 
                                      !!valid the aggregated value.

!used by inverse distance weighting interpolation
REAL (KIND = float), PRIVATE :: idw_power !! power to be used with IDW

!used by kriging
INTEGER (KIND = short), PRIVATE :: krige_var !! when set to 1 a grid of kriging variance is generated and exported if export option is active, default to 0
INTEGER (KIND = short), PRIVATE :: krige_anisotropy !! considers anisotropy, default = 0 excludes anisotropy
INTEGER (KIND = short), PRIVATE :: krige_varmodel !! 1 = spherical, 2 = exponential, 3 = gaussian, 0 = automatic fitting. default to 2
INTEGER (KIND = short), PRIVATE :: krige_lags !! number of lags for semivariogram. if undefined or set to 0 default to 15
REAL (KIND = float)   , PRIVATE :: krige_maxlag !! maximum distance to be considered for semivariogram assessment. If undefined or set to 0, it is computed automatically

INTEGER (KIND = short), PRIVATE :: export !!activates grid exporting
CHARACTER (LEN = 1000), PRIVATE :: export_path  !!folder where to put exported grids
CHARACTER (LEN = 1000), PRIVATE :: export_file  !!name of  exported  grid
CHARACTER (LEN = 1000), PRIVATE :: export_file_var  !!name of  exported  variance grid
TYPE (DateTime), PRIVATE :: export_start   !! time and date to start exporting
TYPE (DateTime), PRIVATE :: export_stop !! time and date to stop exporting
INTEGER (KIND = short), PRIVATE ::  export_dt !! time between two exportations
TYPE (DateTime), PRIVATE :: timeNewExport !! time when new exporting must occur
INTEGER (KIND = short), PRIVATE :: export_format !! 1 = esri_ascii, 2 = esri_binary, 3 = netcdf 
INTEGER (KIND = short), PRIVATE :: export_epsg !! coordinate reference system of exported grid 
TYPE (grid_real), PRIVATE :: exportedGrid, exportedGridVar
TYPE (CRS), PRIVATE :: exportGridMapping
LOGICAL, PRIVATE   :: needConversion



!Public routines
PUBLIC :: AirTemperatureDailyMeanInit
PUBLIC :: AirTemperatureDailyMeanRead


!Local routines
PRIVATE :: SetSpecificProperties


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

!==============================================================================
!| Description:
!   Initialize air temperature
SUBROUTINE AirTemperatureDailyMeanInit &
!
( ini, mask, dtMeteo, tstart, dem_loaded )

IMPLICIT NONE

TYPE (IniList), INTENT(IN) :: ini
TYPE (grid_integer), INTENT(IN) :: mask  !!defines interpolation extent
INTEGER (KIND = short), INTENT(IN) :: dtMeteo !! deltat of meteo data reading
TYPE (DateTime),     INTENT(IN) :: tstart !!initial time
LOGICAL , INTENT (IN) :: dem_loaded !! true if dem has been loaded


!local declarations
CHARACTER (LEN = 1000) :: filename
TYPE (DateTime) :: lapseTime
TYPE (grid_real) :: meteoTemp


!-------------------------end of declarations----------------------------------
!set initial time
  timeNew = tstart  

  !initialize grid
  CALL NewGrid (temperatureAirDailyMean, mask, 0.)
  CALL NewGrid (grid_devst, mask, 0.)
  
  !set deltat
  dtTemperatureDailyMean = day


	!check dt is multiple of dtMeteo
  IF (.NOT.(MOD(dtTemperatureDailyMean,dtMeteo) == 0)) THEN
           CALL Catch ('error', 'AirTemperatureDailyMean',   &
		                 'dt is not multiple of dtMeteo')
  END IF
  
  !set valid threshold
  IF (KeyIsPresent ('valid-threshold', ini, &
                         section = 'temperature-daily-mean') ) THEN
      valid_prcn = IniReadReal ('valid-threshold', ini, &
                                        section = 'temperature-daily-mean')
  ELSE ! set to default = 1.0 
      CALL Catch ('info', 'AirTemperatureDailyMean',   &
		                 'valid-threshold not defined, set to 1')
      valid_prcn = 1.
  END IF

  !set cell-size
  cellsizeInterpolation = mask % cellsize
  
  !interpolation-assignment method
  IF (KeyIsPresent('interpolation-assignment', ini, &
                     section = 'temperature-daily-mean') ) THEN
      interpolationMethod_assignment = IniReadInt ('interpolation-assignment', ini, &
                                    section = 'temperature-daily-mean')
  ELSE
    CALL Catch ('error', 'AirTemperatureDailyMean',   &
                'interpolation-assignment missing in meteo configuration file')
  END IF

  !set interpolation method
  IF (interpolationMethod_assignment == 1) THEN

      interpolationMethod = IniReadInt ('interpolation', ini, &
                                    section = 'temperature-daily-mean')
      
	    CALL SetSpecificProperties (interpolationMethod, ini, mask)
  
  ELSE !read map
      interpolationMethod = -1
      CALL GridByIni (ini, interpolationMethod_map, 'temperature-daily-mean', &
                     "interpolation")
      !scan for interpolation methods included
      interpolationMethod_vector = 0
      DO i = 1, interpolationMethod_map % idim
        DO j = 1, interpolationMethod_map % jdim
            IF (interpolationMethod_map % mat(i,j) /= &
                interpolationMethod_map % nodata) THEN
                interpolationMethod_vector &
                   (interpolationMethod_map % mat (i,j)) = &
                    interpolationMethod_map % mat (i,j)
            END IF
        END DO
      END DO

      DO i = 1, 3
         CALL SetSpecificProperties (interpolationMethod_vector (i), ini, mask)
      END DO

  END IF

  !scale factor and offset
	IF (KeyIsPresent('offset', ini, section = 'temperature-daily-mean')) THEN	
	   offset_value = IniReadReal ('offset', ini, section = 'temperature-daily-mean')
	ELSE
	    offset_value = MISSING_DEF_REAL
	END IF
	
	IF (KeyIsPresent('scale-factor', ini, section = 'temperature-daily-mean')) THEN	
	   scale_factor = IniReadReal ('scale-factor', ini, section = 'temperature-daily-mean')
	ELSE
	    scale_factor = MISSING_DEF_REAL
	END IF			

   !set power_idw
   IF (KeyIsPresent('idw-power', ini, section = 'temperature-daily-mean')) THEN	
	   idw_power = IniReadReal ('idw-power', ini, section = 'temperature-daily-mean')
   ELSE !set default value
	    idw_power = 2.
   END IF	
 
	!file
    filename = IniReadString ('file', ini, section = 'temperature-daily-mean')
    IF (interpolationMethod_assignment == 1 .AND. &
        interpolationMethod == 0) THEN !data are stored in net-cdf file
       !store net-cdf filename
       fileGrid = filename
       IF ( KeyIsPresent ('variable', ini, section = 'temperature-daily-mean') ) THEN
          temperatureAirDailyMean % var_name = IniReadString ('variable', &
                                             ini, section = 'temperature-daily-mean')
            !read grid and store as temporary file
             CALL NewGrid (layer = meteoTemp, fileName = TRIM(fileGrid), &
                         fileFormat = NET_CDF, &
                         variable = TRIM(temperatureAirDailyMean % var_name) )
       ELSE IF  (KeyIsPresent ('standard_name', ini, &
                              section = 'temperature-daily-mean') ) THEN
          temperatureAirDailyMean % standard_name = IniReadString ('standard_name', &
                                              ini, section = 'temperature-daily-mean')
       ELSE
          CALL Catch ('error', 'AirTemperatureDailyMean',   &
		       'standard_name or variable missing in section temeprature' )
       END IF      

       !set cellsize to zero. Optimal cellsize is automatically computed
       cellsizeInterpolation = 0.
            
    ELSE !open file containing local measurements
       fileunit = GetUnit ()
	     OPEN(fileunit,file = filename(1:LEN_TRIM(filename)),status='old')
    END IF
    
  !use elevation to drift interpolation
  IF ( KeyIsPresent ('elevation-drift', ini, section = 'temperature-daily-mean') ) THEN
	     elevationDrift = IniReadInt ('elevation-drift', ini, section = 'temperature-daily-mean')
  ELSE
      elevationDrift = 0 !default, suppress drift
  END IF
  
  IF (elevationDrift == 1) THEN
      
      IF (interpolationMethod == 0 ) THEN
            CALL Catch ('error', 'AirTemperatureDailyMean',   &
                'elevation drift cannot be applied when interpolation = 0 ')
       END IF
      
      !check if have been loaded by domain properties
      IF ( .NOT. dem_loaded) THEN
           CALL Catch ('error', 'AirTemperatureDailyMean',   &
                        'dem for elevation drift was not loaded ')
      END IF
      
       !need lapse-rate computation?
       IF (KeyIsPresent ('lapse-rate-computation', ini, &
                         section = 'temperature-daily-mean') )  THEN
           lapse_rate_computation = IniReadInt ('lapse-rate-computation', ini, &
                                        section = 'temperature-daily-mean')
       ELSE !lapse-rate-computation is missing
           CALL Catch ('error', 'AirTemperatureDailyMean',   &
		                 'lapse-rate-computation missing in meteo configuration file')
       END IF
       
       IF ( lapse_rate_computation == 1 .AND. interpolationMethod == 0 ) THEN                                   
         CALL Catch ('error', 'AirTemperatureDailyMean',   &
            'lapse_rate_computation option requires interpolation different from 0')
       END IF
       
       IF (KeyIsPresent ('r2min', ini, &
                         section = 'temperature-daily-mean') )  THEN
           r2min = IniReadInt ('r2min', ini, &
                                        section = 'temperature-daily-mean')
       ELSE !r2min is missing
            IF (lapse_rate_computation == 1) THEN                                   
                 CALL Catch ('error', 'AirTemperatureDailyMean',   &
		                 'r2min missing in meteo configuration file')
            END IF
       END IF
      
       
        !set lapse rate assignment
        IF (KeyIsPresent ('lapse-rate-assignment', ini, &
                            section = 'temperature-daily-mean') )  THEN
            lapse_rate_assignment = IniReadInt ('lapse-rate-assignment', ini, &
                                        section = 'temperature-daily-mean')
        ELSE !lapse-rate-assignment missing
            IF (lapse_rate_computation == 0) THEN         
                CALL Catch ('error', 'AirTemperatureDailyMean',   &
                'lapse-rate-assignment missing in meteo configuration file')
            ELSE
                lapse_rate_assignment = 1
            END IF
        END IF
       
        !lapse rate map
        IF (lapse_rate_assignment == 2) THEN
            IF ( SubSectionIsPresent (subsection = 'lapse-rate-map', &
                    section = 'temperature-daily-mean', iniDb = ini) ) THEN

                !check all keywords are defined
                IF (IniReadString ('format', ini, 'temperature-daily-mean',&
                                        'lapse-rate-map') == 'net-cdf') THEN

                    IF (KeyIsPresent ('time',ini, 'temperature-daily-mean', &
                                                    'lapse-rate-map'))  THEN
                        timeString = IniReadString ('time', ini, &
                                'temperature-daily-mean', 'lapse-rate-map') 
                        lapseTime = timeString 
                        IF (lapseTime > tstart) THEN
                            CALL Catch ('error', 'AirTemperatureDailyMean',   &
		                    'lapse-rate-map time greater than starting time' )
                        END IF
                    ELSE
                        CALL Catch ('error', 'AirTemperatureDailyMean',   &
		                    'time keyword is missing in lapse-rate-map subsection' )
                    END IF

                END IF
              
                CALL GridByIni (ini, lapse_map, 'temperature-daily-mean', "lapse-rate-map") 
                !check coordinate reference system and spatial resolution
                IF (.NOT. CRSisEqual (mask = mask, grid = lapse_map, &
                                    checkCells= .TRUE.)) THEN
                    CALL Catch ('error', 'AirTemperatureDailyMean',   &
		                    'wrong lapse-rate-map spatial reference system ')
                END IF
            ELSE
                CALL Catch ('error', 'AirTemperatureDailyMean',   &
		                    'lapse-rate-map missing in meteo configuration file')
            END IF
        ELSE !build lapse map from scalar value

            IF ( KeyIsPresent ('lapse-rate-scalar', ini, &
                    section = 'temperature-daily-mean') ) THEN
                lapse = IniReadReal ('lapse-rate-scalar', ini, &
                    section = 'temperature-daily-mean')
                CALL NewGrid (lapse_map, mask)
                lapse_map = lapse
            ELSE
                CALL Catch ('error', 'AirTemperatureDailyMean',   &
		                    'lapse-rate-scalar missing in meteo configuration file')
            END IF
        END IF
		   
	END IF  



  !grid exporting settings
  IF (KeyIsPresent ('export', ini, section = 'temperature-daily-mean')  )  THEN
     export = IniReadInt ('export', ini, section = 'temperature-daily-mean')
  ELSE
     export = 0
  END IF 

  IF (export == 1) THEN
     
     
     !set export path name
     IF (KeyIsPresent ('export-path', ini, section = 'temperature-daily-mean')  )  THEN
         export_path = IniReadString ('export-path', ini, section = 'temperature-daily-mean')
         !check if path ends with / or \
         IF (  export_path (LEN_TRIM (export_path):LEN_TRIM (export_path)) /= '\' .AND. &
               export_path (LEN_TRIM (export_path):LEN_TRIM (export_path)) /= '/' ) THEN
             export_path (LEN_TRIM (export_path)+1:LEN_TRIM (export_path)+1) = '\'
         END IF
       
         IF (INDEX (export_path,'.') == 1) THEN !detected relative path, remove '.'
            export_path = export_path (2:LEN_TRIM(export_path))
            !build absolute path
            export_path = TRIM(CurrentDir() ) // TRIM(export_path)
         END IF

         !check OS convention
        IF (GetOS () == WIN32) THEN
          export_path = ReplaceChar (export_path,'/','\')
        ELSE
          export_path = ReplaceChar (export_path,'\','/')
        END IF
  
         
         !check folder exists
         IF ( .NOT. DirExists (TRIM (export_path) ) ) THEN
             CALL Catch ('info', 'AirTemperatureDailyMean',   &
                  'creating directory:  ', argument = TRIM(export_path))
             CALL DirNew (export_path)
         END IF
     ELSE
         CALL Catch ('error', 'AirTemperatureDailyMean',   &
                  'missing export-path')
     END IF 
     
     !starting time
     IF (KeyIsPresent ('export-start', ini, section = 'temperature-daily-mean')  )  THEN
         timeString = IniReadString ('export-start', ini, 'temperature-daily-mean')
         export_start = timeString
     ELSE
         CALL Catch ('error', 'AirTemperatureDailyMean',   &
                  'missing export-start')
     END IF 
     
     !initialize timeNewExport
     timeNewExport = export_start
     
     !ending time
     IF (KeyIsPresent ('export-stop', ini, section = 'temperature-daily-mean')  )  THEN
         timeString = IniReadString ('export-stop', ini, 'temperature-daily-mean')
         export_stop = timeString
     ELSE
         CALL Catch ('error', 'AirTemperatureDailyMean',   &
                  'missing export-start')
     END IF 
     
     !export dt
     IF (KeyIsPresent ('export-dt', ini, section = 'temperature-daily-mean')  )  THEN
         export_dt = IniReadInt ('export-dt', ini, section = 'temperature-daily-mean')
     ELSE
         CALL Catch ('error', 'AirTemperatureDailyMean',   &
                  'missing export-dt')
     END IF 

     !coordinate reference system
     IF (KeyIsPresent ('export-epsg', ini, section = 'temperature-daily-mean')  )  THEN
         export_epsg = IniReadInt ('export-epsg', ini, section = 'temperature-daily-mean')
         
         exportGridMapping = DecodeEPSG (export_epsg)
         IF (exportGridMapping == temperatureAirDailyMean % grid_mapping) THEN
            needConversion = .FALSE.
            !initialize grid for exporting with CRS as meteo variable
            CALL NewGrid (exportedGrid, temperatureAirDailyMean)
            CALL NewGrid (exportedGridVar, temperatureAirDailyMean)
         ELSE
            needConversion = .TRUE.
            !initialize grid for converting with required CRS
            exportedGrid % grid_mapping = DecodeEPSG (export_epsg)
            exportedGridVar % grid_mapping = DecodeEPSG (export_epsg)
         END IF
     ELSE
         CALL Catch ('info', 'AirTemperatureDailyMean',   &
                  'export-epsg not defined, use default')
         needConversion = .FALSE.
         !initialize grid for exporting with CRS 
         CALL NewGrid (exportedGrid, temperatureAirDailyMean)
         CALL NewGrid (exportedGridVar, temperatureAirDailyMean)
     END IF 

     !export file format 
     IF (KeyIsPresent ('export-format', ini, section = 'temperature-daily-mean')  )  THEN
         export_format = IniReadInt ('export-format', ini, section = 'temperature-daily-mean')
     ELSE
         CALL Catch ('error', 'AirTemperatureDailyMean',   &
                  'missing export-format')
     END IF 

     IF (export_format == 3) THEN
       !grid map  
       CALL SetLongName ( 'air_temperature', exportedGrid)
       CALL SetStandardName ( 'air_temperature', exportedGrid)
       CALL SetUnits ('degree_Celsius', exportedGrid) !this unit is for exporting grid, it is converted internally to m/s
       !if file exists, remove it
       export_file = TRIM(export_path) //  'temperature_daily_mean.nc'
       IF ( FileExists (export_file) ) THEN
          CALL FileDelete (export_file)
       END IF
       
       !variance of kriging 
       CALL SetLongName ( 'air_temperature_variance', exportedGridVar)
       CALL SetStandardName ( 'air_temperature_variance', exportedGridVar)
       CALL SetUnits ('degree_Celsius', exportedGridVar) !this unit is for exporting grid
       !if file exists, remove it
       export_file_var = TRIM(export_path) //  'temperature_daily_mean_variance.nc'
       IF ( FileExists (export_file_var) ) THEN
          CALL FileDelete (export_file_var)
       END IF
       
     END IF

  END IF
	
	!complete initialization

  IF (interpolationMethod == 0) THEN
				!Get the dt of imported  field. Assume dt is regular	
				dtGrid = GetDtGrid (filename = fileGrid, checkRegular = .TRUE.)
				!check dt is multiple of dtGrid
        IF (.NOT.(MOD(dtTemperatureDailyMean,dtGrid) == 0)) THEN
            CALL Catch ('error', 'AirTemperatureDailyMean',   &
                'dt temperature daily mean is not multiple of dt of input grid')
        END IF
   ELSE
        !populate  metadata
        CALL ReadMetadata (fileunit, thermometersDailyMean)

        !check spatial reference system
        IF ( .NOT. thermometersDailyMean % mapping == mask % grid_mapping)  THEN
            CALL Catch ('info', 'AirTemperatureDailyMean',   &
		              'converting coordinate of stations')
            !convert stations' coordinate
            point1 % system = thermometersDailyMean % mapping
            point2 % system = mask % grid_mapping
            thermometersDailyMean % mapping = mask % grid_mapping
            DO i = 1, thermometersDailyMean % countObs
              point1 % easting = thermometersDailyMean % obs (i) % xyz % easting
              point1 % northing = thermometersDailyMean % obs (i) % xyz % northing
              point1 % elevation = thermometersDailyMean % obs (i) % z
              CALL Convert (point1, point2)
              thermometersDailyMean % obs (i) % xyz % easting = point2 % easting
              thermometersDailyMean % obs (i) % xyz % northing = point2 % northing 
            END DO
        END IF
                         
        !set supplementary stations network for elevation drift
        IF (elevationDrift == 1) THEN
          stationsRefElev =  thermometersDailyMean
          DO i = 1, thermometersDailyMean % countObs
              stationsRefElev % obs (i) % xyz % elevation = refElevation
          END DO     
	
        END IF
        
        !allocate vectors if lapse rate has to be computed from data
        IF ( lapse_rate_computation == 1) THEN
            ALLOCATE (vectorT (thermometersDailyMean % countObs) )
            ALLOCATE (vectorZ (thermometersDailyMean % countObs))
        END IF
             
  END IF
   
RETURN
END SUBROUTINE AirTemperatureDailymeanInit

!==============================================================================
!| Description:
!   Read air temperature data
SUBROUTINE AirTemperatureDailyMeanRead &
!
( time, dem )

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT(IN) :: time !!current time
TYPE (grid_real), INTENT(IN) :: dem  !digital elevation model to apply drift

!local declarations:
TYPE (DateTime) :: time_toread !! time to start reading from
CHARACTER (LEN = 300) :: filename
CHARACTER (LEN = 300) :: varname
REAL (KIND = float)   :: rsquare
INTEGER (KIND = short) :: i, j

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

IF ( .NOT. (time < timeNew) ) THEN
   
   !time_toread = time + - (dtTemperatureDailyMean - thermometersDailyMean % timeIncrement)
   time_toread = time
   timeString = time_toread
   CALL Catch ('info', 'AirTemperatureDailyMean',   &
		                 'read new temperature data: ', &
                     argument = timeString)
   
   !update lapse rate from maps when required
    IF (elevationDrift == 1) THEN
       IF (time == lapse_map % next_time) THEN !update lapse_map
          varname = lapse_map %var_name
          filename = lapse_map % file_name
          !destroy current grid
          CALL GridDestroy (lapse_map)
          !read new grid in netcdf file
          CALL NewGrid (layer = lapse_map, fileName = TRIM(filename), &
                         fileFormat = NET_CDF, &
                         variable = TRIM(varname), &
                         time = time_toread)
      END IF
    END IF

   SELECT CASE (interpolationMethod)
    CASE (0) !input grid
      
        CALL ReadField (fileGrid,  time_toread, &
                        dtTemperatureDailyMean, dtGrid, &
                        'M', temperatureAirDailyMean, &
                         varName = temperatureAirDailyMean % var_name)
 
    CASE DEFAULT !requires interpolation
      !read new  data
      CALL ReadData (network = thermometersDailyMean, fileunit = fileunit, &
                      time = time_toread, aggr_time = dtTemperatureDailyMean, &
                      aggr_type = 'ave', tresh = valid_prcn)
      
      IF (elevationDrift == 1) THEN
        !update lapse rate from data  
        IF (lapse_rate_computation == 1) THEN
            
          DO i = 1, thermometersDailyMean % countObs
               vectorT (i) = thermometersDailyMean % obs (i) % value
               vectorZ (i) = thermometersDailyMean % obs (i) % xyz % elevation
          END DO
          
          lapse_computed = LinearRegressionSlope (x=vectorZ, y=vectorT, &
                    nodata = thermometersDailyMean % nodata, r2 = rsquare)
          IF (rsquare < r2min) THEN
              lapse_computed = lapse
          ENDIF
          
          lapse_map = lapse_computed
        END IF

        !shift  observation to reference elevation by applying lapse rate
        CALL ShiftMeteoWithLapse (thermometersDailyMean, lapse_map, refElevation, &
                                  stationsRefElev)
        
        !interpolate 
        IF (interpolationMethod_assignment == 1) THEN !only one method is applied
                
                CALL Interpolate (network = stationsRefElev, &
                            method = interpolationMethod, &
                            near = neighbors, &
                            idw_power = idw_power, & 
                            anisotropy = krige_anisotropy, &
                            varmodel = krige_varmodel, &
                            lags = krige_lags, &
                            maxlag = krige_maxlag, &
                            grid = temperatureAirDailyMean, &
                            devst = grid_devst)                  
                
        ELSE 
            !loop trough interpolation methods
            DO i = 1, 3
                IF (interpolationMethod_vector(i) > 0) THEN
                      
                    CALL Interpolate (network = stationsRefElev, &
                            method = interpolationMethod_vector(i), &
                            near = neighbors, &
                            idw_power = idw_power, & 
                            anisotropy = krige_anisotropy, &
                            varmodel = krige_varmodel, &
                            lags = krige_lags, &
                            maxlag = krige_maxlag, &
                            grid = interpolatedMap (i), &
                            devst = grid_devst) 
                      
                END IF
            END DO

            !compose the final interpolated field
            DO i = 1, interpolationMethod_map % idim
                DO j = 1, interpolationMethod_map % jdim
                    IF (interpolationMethod_map % mat(i,j) /= &
                        interpolationMethod_map % nodata ) THEN
                        temperatureAirDailyMean % mat (i,j) = &
                        interpolatedMap (interpolationMethod_map % mat(i,j)) % mat(i,j)
                    END IF
                END DO
            END DO
        END IF

        !Shift back interpolated field to terrain surface
        CALL ShiftBackWithLapse (temperatureAirDailyMean, dem, &
                                lapse_map, refElevation)
       
         
      ELSE !elevationdrift = 0
        

        IF (interpolationMethod_assignment == 1) THEN !only one method is applied
                
                CALL Interpolate (network = thermometersDailyMean, &
                            method = interpolationMethod, &
                            near = neighbors, &
                            idw_power = idw_power, & 
                            anisotropy = krige_anisotropy, &
                            varmodel = krige_varmodel, &
                            lags = krige_lags, &
                            maxlag = krige_maxlag, &
                            grid = temperatureAirDailyMean, &
                            devst = grid_devst) 
                
        ELSE
            !loop trough interpolation methods
            DO i = 1, 3
                IF (interpolationMethod_vector(i) > 0) THEN
                      
                    CALL Interpolate (network = thermometersDailyMean, &
                            method = interpolationMethod_vector(i), &
                            near = neighbors, &
                            idw_power = idw_power, & 
                            anisotropy = krige_anisotropy, &
                            varmodel = krige_varmodel, &
                            lags = krige_lags, &
                            maxlag = krige_maxlag, &
                            grid =  interpolatedMap (i), &
                            devst = grid_devst) 
                    
                END IF
            END DO

            !compose the final interpolated field
            DO i = 1, interpolationMethod_map % idim
                DO j = 1, interpolationMethod_map % jdim
                    IF (interpolationMethod_map % mat(i,j) /= &
                        interpolationMethod_map % nodata ) THEN
                        temperatureAirDailyMean % mat (i,j) = &
                        interpolatedMap (interpolationMethod_map % mat(i,j)) % mat(i,j)
                    END IF
                END DO
            END DO

        END IF  !1 or more interpolation methods       
      END IF  !elevationDrift
		
	END SELECT

  !apply scale factor and offset, if defined
	IF (offset_value /= MISSING_DEF_REAL) THEN
	   CALL Offset (temperatureAirDailyMean, offset_value)
	END IF
	
	
	IF (scale_factor /= MISSING_DEF_REAL) THEN
	   CALL Scale (temperatureAirDailyMean, scale_factor)
	END IF


  !grid exporting
  IF(export > 0 ) THEN
	  IF( time == timeNewExport .AND. &
        time >= export_start .AND. &
        time <= export_stop ) THEN
        timeString = time
        timeString = timeString(1:19)
        timeString(14:14) = '-'
		    timeString(17:17) = '-'
        
        grid_devst % reference_time = temperatureAirDailyMean % reference_time
        IF (needConversion) THEN
           CALL GridConvert (temperatureAirDailyMean, exportedGrid)
           CALL GridConvert (grid_devst, exportedGridVar)
        ELSE
           exportedGrid = temperatureAirDailyMean
           exportedGridVar = grid_devst
        END IF 

        SELECT CASE (export_format)
        CASE (1) !esri-ascii
              export_file = TRIM(export_path) //  TRIM(timeString) // &
                           '_temperature_daily_mean.asc'
              CALL Catch ('info', 'AirTemperatureDailyMean',   &
		                       'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGrid, export_file, ESRI_ASCII)
              
              IF (krige_var == 1) THEN !export kriging variance
                    export_file_var = TRIM(export_path) //  TRIM(timeString) // &
                           '_temperature_daily_mean_variance.asc'
                    CALL Catch ('info', 'AirTemperatureDailyMean',   &
		                       'exporting variance grid to file: ', &
                              argument = TRIM(export_file_var))
                    CALL ExportGrid (exportedGridVar, export_file_var, ESRI_ASCII)
              END IF
              
        CASE (2) !esri-binary
              export_file = TRIM(export_path) //  TRIM(timeString) // &
                           '_temperature_daily_mean'
              CALL Catch ('info', 'AirTemperatureDailyMean',   &
		                       'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGrid, export_file, ESRI_BINARY)
              
              IF (krige_var == 1) THEN !export kriging variance
                   export_file_var = TRIM(export_path) //  TRIM(timeString) // &
                               '_temperature_daily_mean_variance'
                   CALL Catch ('info', 'AirTemperatureDailyMean',   &
		                       'exporting variance grid to file: ', &
                              argument = TRIM(export_file_var))
                   CALL ExportGrid (exportedGridVar, export_file_var, ESRI_BINARY)
              END IF
              
        CASE (3) !net_cdf
              CALL SetCurrentTime (time, exportedGrid)
              CALL Catch ('info', 'AirTemperatureDailyMean',   &
                           'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGrid, export_file, 'air_temperature', 'append')
              
              IF (krige_var == 1) THEN !export kriging variance
                  CALL SetCurrentTime (time, exportedGridVar)
                  CALL Catch ('info', 'AirTemperatureDailyMean',   &
		                       'exporting grid to file: ', &
                              argument = TRIM(export_file_var))
                 CALL ExportGrid (exportedGridVar, export_file_var, 'air_temperature_variance', 'append')
              END IF
        END SELECT
        timeNewExport = timeNewExport + export_dt
    END IF
  ENDIF



  !time forward
  timeNew = timeNew + dtTemperatureDailyMean
END IF

RETURN
END SUBROUTINE AirTemperatureDailyMeanRead


!==============================================================================
!! Description:
!!   set properties and initialize variables for each interpolation method
SUBROUTINE SetSpecificProperties &
!
( method, ini, mask )

IMPLICIT NONE

!Arguments with intent(in):
INTEGER (KIND = short), INTENT(in) :: method
TYPE (IniList)        , INTENT(in) :: ini
TYPE (grid_integer)   , INTENT(in) :: mask

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

!set data specific for each interpolation method
	SELECT CASE (method)
    CASE (1) !thiessen
      CALL NewGrid (interpolatedMap (1), mask, 0.)
    CASE (2) !inverse distance
      CALL NewGrid (interpolatedMap (2), mask, 0.)
	    !nearest-points
      IF (KeyIsPresent ('nearest-points', ini, &
                         section = 'temperature-daily-mean') ) THEN
	       neighbors = IniReadInt ('nearest-points', ini, &
                                   section = 'temperature-daily-mean')
      ELSE !nearest-points missing
          CALL Catch ('error', 'AirTemperatureDailyMean',   &
		                 'nearest-points missing in meteo configuration file')
      END IF
 
    CASE (3) !kriging
      CALL NewGrid (interpolatedMap (3), mask, 0.)
      !nearest-points
      IF (KeyIsPresent ('nearest-points', ini, &
                         section = 'temperature-daily-mean') ) THEN
	       neighbors = IniReadInt ('nearest-points', ini, &
                                   section = 'temperature-daily-mean')
           CALL Catch ('info', 'AirTemperatureDailyMean', &
               'neighbors set to: ', argument = ToString (neighbors) )
      ELSE !nearest-points missing
          CALL Catch ('error', 'AirTemperatureDailyMean',   &
		                 'nearest-points missing in meteo configuration file')
      END IF
      
      !kriging-variance
      IF (KeyIsPresent ('kriging-variance', ini, &
                         section = 'temperature-daily-mean') ) THEN
	       krige_var = IniReadInt ('kriging-variance', ini, &
                                   section = 'temperature-daily-mean')
           
           CALL Catch ('info', 'AirTemperatureDailyMean',   &
		                 'kriging variance export set to ', argument = ToString(krige_var) )
      ELSE !kriging-variance missing
           krige_var = 0                                     
           CALL Catch ('info', 'AirTemperatureDailyMean',   &
		                 'kriging variance export set to default (0)')
      END IF
      
      !kriging-anisotropy
      IF (KeyIsPresent ('kriging-anisotropy', ini, &
                         section = 'temperature-daily-mean') ) THEN
	       krige_anisotropy = IniReadInt ('kriging-anisotropy', ini, &
                                   section = 'temperature-daily-mean')
           CALL Catch ('info', 'AirTemperatureDailyMean', &
               'kriging anisotropy set to: ', argument = ToString (krige_anisotropy) )
      ELSE !kriging-anisotropy missing
           krige_anisotropy = 0                                     
           CALL Catch ('info', 'AirTemperatureDailyMean',   &
		                 'kriging anisotropy set to default (0)')
      END IF
      
      !kriging-variogram-model
      IF (KeyIsPresent ('kriging-variogram-model', ini, &
                         section = 'temperature-daily-mean') ) THEN
	       krige_varmodel = IniReadInt ('kriging-variogram-model', ini, &
                                   section = 'temperature-daily-mean')
           CALL Catch ('info', 'AirTemperatureDailyMean', &
               'kriging variogram model set to: ', argument = ToString (krige_varmodel) )
           IF (krige_varmodel == 0 ) krige_varmodel = 5 !automatic fitting
      ELSE !kriging-variogram-model
           krige_varmodel = 2                                     
           CALL Catch ('info', 'AirTemperatureDailyMean',   &
		                 'kriging variogram model set to default (2 exponential)')
      END IF
      
      !kriging-lags
      IF (KeyIsPresent ('kriging-lags', ini, &
                         section = 'temperature-daily-mean') ) THEN
	       krige_lags = IniReadInt ('kriging-lags', ini, &
                                   section = 'temperature-daily-mean')
           IF (krige_lags == 0) krige_lags = 15
           CALL Catch ('info', 'AirTemperatureDailyMean', &
               'kriging lags set to: ', argument = ToString (krige_lags) )
      ELSE !kriging-variogram-model
           krige_lags = 15                                     
           CALL Catch ('info', 'AirTemperatureDailyMean',   &
		                 'kriging lags set to default (15)')
      END IF
      
       !kriging-maxlag
      IF (KeyIsPresent ('kriging-maxlag', ini, &
                         section = 'temperature-daily-mean') ) THEN
	       krige_maxlag = IniReadInt ('kriging-maxlag', ini, &
                                   section = 'temperature-daily-mean')
           CALL Catch ('info', 'AirTemperatureDailyMean', &
               'kriging max lag set to: ', argument = ToString (krige_maxlag) )
      ELSE !kriging-maxlag missing
           krige_maxlag = 0
          CALL Catch ('error', 'AirTemperatureDailyMean',   &
		                 'kriging max lag set to default (0)')
      END IF
	
	END SELECT


RETURN
END SUBROUTINE SetSpecificProperties


END MODULE AirTemperatureDailyMean