!! Manage wind speed data with arbitrary time cumulation
!|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 - 4th June 2021  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 4/Jun/2021 | Original code |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! Module to manage wind speed data with arbitrary time cumulation
!
MODULE WindFlux


! 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 Morphology, ONLY : &
    !Imported routines:
    CurvatureMicroMet, DeriveSlope, &
    DeriveAspect

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, Interpolate, &
  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 SpatialInterpolation, ONLY: &
    !Imported routines:
    WindMicromet, WindGonzalezLongatt
    
IMPLICIT NONE 

! Global (i.e. public) Declarations: 

INTEGER (KIND = short)     :: dtWindSpeed  = 0!!cumulation deltat 
TYPE(ObservationalNetwork) :: anemometersWS !!wind speed stations network
TYPE(grid_real)            :: windSpeed	   	  !![m/s]


!Private declarations
TYPE(DateTime), PRIVATE         :: timeNew !!time when new  data must be read
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 (5) !!defines active interpolation methods 
TYPE (grid_real),       PRIVATE :: interpolatedMap (5)  !!1 = map for thiessen, 2 = map for idw, &
                                                        !3 = map for kriging, 4 = map for micromet, &
                                                        !5 = map for gonzalez
REAL (KIND = float),    PRIVATE :: cellsizeInterpolation !!spatial resolution of interpolated grid
INTEGER (KIND = short), PRIVATE :: neighbors  !!number of closest data to use for interpolation
TYPE(grid_real),        PRIVATE :: grid_devst !!standard deviation of kriging interpolation
INTEGER (KIND = short), PRIVATE :: elevationDrift !! 1 = use elevation to modify interpolated data
LOGICAL,                PRIVATE :: micrometInUse !! algorithm micromet in use for interpolation
INTEGER (KIND = short), PRIVATE :: timezone !!local timezone
TYPE (grid_real),       PRIVATE :: slope !terrain slope grid, used by micromet
TYPE (grid_real),       PRIVATE :: curvature !terrain curvature grid, used by micromet
TYPE (grid_real),       PRIVATE :: slope_az !terrain slope azimuth (aspect), used by micromet
REAL (KIND = float),    PRIVATE :: micrometSlopeWT !! Micromet slope weighting factor. default = 0.5
REAL (KIND = float),    PRIVATE :: micrometCurvatureWT !! Micromet curvature weighting factor. default = 0.5
REAL (KIND = float),    PRIVATE :: micrometLengthScale !! length scale used by micromet to compute curvature. default = 5000 m
TYPE(ObservationalNetwork), PRIVATE :: anemometersWD !!wind direction stations network
INTEGER (KIND = short), PRIVATE :: fileunitWD !! unit of file containing wind direction data
CHARACTER (LEN = 300),  PRIVATE :: fileGrid !!file containing grids
INTEGER (KIND = short), PRIVATE :: dtGrid !!dt of imported  field
INTEGER (KIND = short), PRIVATE :: i, j  !!loop index
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
!INTEGER (KIND = short), PRIVATE :: exportGridPrecipitation
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 :: WindFluxInit
PUBLIC :: WindFluxRead


!Local routines
PRIVATE :: SetSpecificProperties


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

!==============================================================================
!| Description:
!   Initialize air temperature
SUBROUTINE WindFluxInit &
!
( ini, mask, dtMeteo, tstart, dem, 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
TYPE(grid_real),    INTENT(in)  :: dem !!digital elevation model to be used to modify interpolated data
LOGICAL , INTENT (IN) :: dem_loaded !! true if dem has been loaded


!local declarations
CHARACTER (LEN = 1000) :: filename, filenameWD
TYPE (grid_real) :: meteoTemp


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

  !initialize grid
  CALL NewGrid (windSpeed, mask, 0.)
  CALL NewGrid (grid_devst, mask, 0.)
  
  !set deltat
  dtWindSpeed = IniReadInt ('dt', ini, section = 'wind-speed')

  !check dt is multiple of dtMeteo
  IF (.NOT.(MOD(dtWindSpeed,dtMeteo) == 0)) THEN
           CALL Catch ('error', 'WindFlux',   &
		                 'dt is not multiple of dtMeteo')
  END IF
  
  !set valid threshold
  IF (KeyIsPresent ('valid-threshold', ini, &
                         section = 'wind-speed') ) THEN
      valid_prcn = IniReadReal ('valid-threshold', ini, &
                                        section = 'wind-speed')
  ELSE ! set to default = 1.0 
      CALL Catch ('info', 'WindFlux',   &
		                 '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 = 'wind-speed') ) THEN
      interpolationMethod_assignment = IniReadInt ('interpolation-assignment', ini, &
                                    section = 'wind-speed')
  ELSE
    CALL Catch ('error', 'WindFlux',   &
                'interpolation-assignment missing in meteo configuration file')
  END IF   

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

      interpolationMethod = IniReadInt ('interpolation', ini, &
                                    section = 'wind-speed')
      
	    CALL SetSpecificProperties (interpolationMethod, ini, mask)
  
  ELSE !read map
      interpolationMethod = -1
      CALL GridByIni (ini, interpolationMethod_map, 'wind-speed', &
                     '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, 5
         CALL SetSpecificProperties (interpolationMethod_vector (i), ini, mask)
      END DO

  END IF

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

   !set power_idw
   IF (KeyIsPresent('idw-power', ini, section = 'wind-speed')) THEN	
	   idw_power = IniReadReal ('idw-power', ini, section = 'wind-speed')
   ELSE !set default value
	    idw_power = 2.
   END IF	
 
	!file
    filename = IniReadString ('file', ini, section = 'wind-speed')
    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 = 'wind-speed') ) THEN
          windSpeed % var_name = IniReadString ('variable', &
                                             ini, section = 'wind-speed')
            !read grid and store as temporary file
             CALL NewGrid (layer = meteoTemp, fileName = TRIM(fileGrid), &
                         fileFormat = NET_CDF, &
                         variable = TRIM(windSpeed % var_name) )
       ELSE IF  (KeyIsPresent ('standard_name', ini, &
                              section = 'wind-speed') ) THEN
          windSpeed % standard_name = IniReadString ('standard_name', &
                                              ini, section = 'wind-speed')
       ELSE
          CALL Catch ('error', 'WindFlux',   &
		       'standard_name or variable missing in section wind-speed' )
       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
    
    elevationDrift = 0
    micrometInUse = .FALSE.
    !set elevationfridt option
    IF ( interpolationMethod_assignment == 1 ) THEN
        
        IF ( interpolationMethod == 4 ) THEN  !micromet
           elevationDrift = 1
           micrometInUse = .TRUE.
        END IF
        
        IF ( interpolationMethod == 5 ) THEN  !gonzalez
            elevationDrift = 1
        END IF
            
    ELSE !interpolationMethod_assignment = 2
        !parse interpolationMethod_map
        DO i = 1, interpolationMethod_map % idim
            DO j = 1, interpolationMethod_map % jdim
                IF ( interpolationMethod_map % mat (i,j) == 4 ) THEN ! micromet
                     elevationDrift = 1
                     micrometInUse = .TRUE.
                END IF
                IF ( interpolationMethod_map % mat (i,j) == 5 ) THEN !gonzalez
                    elevationDrift = 1
                END IF
            END DO
        END DO 
    END IF
      
    IF (elevationDrift == 1 ) THEN
      
          !check if dem have been loaded by domain properties
          IF ( .NOT. dem_loaded) THEN
               CALL Catch ('error', 'WindFlux',   &
                            'dem for elevation drift was not loaded ')
          END IF
          
          
          !load wind direction metadata
         IF (KeyIsPresent ('wind-direction-file', ini, &
                         section = 'wind-speed') ) THEN
               filenameWD = IniReadString ('wind-direction-file', &
                                       ini, section = 'wind-speed')
        ELSE
             CALL Catch ('error', 'WindFlux',   &
		                 'missing wind speed direction file &
                          in section wind-speed')                
        
        END IF
          
        fileunitWD = GetUnit ()
        
        OPEN (fileunitWD, file = filenameWD(1:LEN_TRIM(filenameWD)), status='old')

        CALL ReadMetadata (fileunitWD, anemometersWD)
          
    END IF
    
    IF ( micrometInUse ) THEN
       
       !load micromet parameters values
       IF (KeyIsPresent ('micromet-length-scale', ini, &
                         section = 'wind-speed') ) THEN
           micrometLengthScale = IniReadReal ('micromet-length-scale', ini, &
                                        section = 'wind-speed')
       ELSE ! set to default = 5000
            CALL Catch ('info', 'WindFlux',   &
		                 'micrometLengthScale set to 5000 m')
            micrometLengthScale = 5000.
       END IF
       
       IF (KeyIsPresent ('micromet-slopewt', ini, &
                         section = 'wind-speed') ) THEN
           micrometSlopeWT = IniReadReal ('micromet-slopewt', ini, &
                                        section = 'wind-speed')
       ELSE ! set to default = 0.5
            CALL Catch ('info', 'WindFlux',   &
		                 'micrometSlopeWT set to 0.5')
            micrometSlopeWT = 0.5
       END IF
       
       IF (KeyIsPresent ('micromet-curvewt', ini, &
                         section = 'wind-speed') ) THEN
           micrometCurvatureWT = IniReadReal ('micromet-curvewt', ini, &
                                        section = 'wind-speed')
       ELSE ! set to default = 0.5
            CALL Catch ('info', 'WindFlux',   &
		                 'micrometCurvatureWT set to 0.5')
            micrometCurvatureWT = 0.5
       END IF
       
       !Compute the curvature
       CALL CurvatureMicroMet (dem, micrometLengthScale, curvature)

       !Compute the terrain slope [rad]
       CALL DeriveSlope (dem, slope)
  
       ! Compute the terrain slope azimuth [rad]
       CALL DeriveAspect (dem, slope_az)
       
        
    END IF  

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

  IF (export == 1) THEN
     
     
     !set export path name
     IF (KeyIsPresent ('export-path', ini, section = 'wind-speed')  )  THEN
         export_path = IniReadString ('export-path', ini, section = 'wind-speed')
         !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', 'WindFlux',   &
                  'creating directory:  ', argument = TRIM(export_path))
             CALL DirNew (export_path)
         END IF
     ELSE
         CALL Catch ('error', 'WindFlux',   &
                  'missing export-path')
     END IF 
     
     !starting time
     IF (KeyIsPresent ('export-start', ini, section = 'wind-speed')  )  THEN
         timeString = IniReadString ('export-start', ini, 'wind-speed')
         export_start = timeString
     ELSE
         CALL Catch ('error', 'WindFlux',   &
                  'missing export-start')
     END IF 
     
     !initialize timeNewExport
     timeNewExport = export_start
     
     !ending time
     IF (KeyIsPresent ('export-stop', ini, section = 'wind-speed')  )  THEN
         timeString = IniReadString ('export-stop', ini, 'wind-speed')
         export_stop = timeString
     ELSE
         CALL Catch ('error', 'WindFlux',   &
                  'missing export-start')
     END IF 
     
     !export dt
     IF (KeyIsPresent ('export-dt', ini, section = 'wind-speed')  )  THEN
         export_dt = IniReadInt ('export-dt', ini, section = 'wind-speed')
     ELSE
         CALL Catch ('error', 'WindFlux',   &
                  'missing export-dt')
     END IF 

     !coordinate reference system
     IF (KeyIsPresent ('export-epsg', ini, section = 'wind-speed')  )  THEN
         export_epsg = IniReadInt ('export-epsg', ini, section = 'wind-speed')
         
         exportGridMapping = DecodeEPSG (export_epsg)
         IF (exportGridMapping == windSpeed % grid_mapping) THEN
            needConversion = .FALSE.
            !initialize grid for exporting with CRS as meteo variable
            CALL NewGrid (exportedGrid, windSpeed)
            CALL NewGrid (exportedGridVar, windSpeed)
         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', 'WindFlux',   &
                  'export-epsg not defined, use default')
         needConversion = .FALSE.
         !initialize grid for exporting with CRS 
         CALL NewGrid (exportedGrid, windSpeed)
         CALL NewGrid (exportedGridVar, windSpeed)
     END IF 

     !export file format 
     IF (KeyIsPresent ('export-format', ini, section = 'wind-speed')  )  THEN
         export_format = IniReadInt ('export-format', ini, section = 'wind-speed')
     ELSE
         CALL Catch ('error', 'WindFlux',   &
                  'missing export-format')
     END IF 

     IF (export_format == 3) THEN
       !grid map  
       CALL SetLongName ( 'wind_speed', exportedGrid)
       CALL SetStandardName ( 'wind_speed', exportedGrid)
       CALL SetUnits ('m/s', exportedGrid) 
       !if file exists, remove it
       export_file = TRIM(export_path) //  'wind_speed.nc'
       IF ( FileExists (export_file) ) THEN
          CALL FileDelete (export_file)
       END IF
       
       !variance of kriging 
       CALL SetLongName ( 'wind_speed_variance', exportedGridVar)
       CALL SetStandardName ( 'wind_speed_variance', exportedGridVar)
       CALL SetUnits ('m/s', exportedGridVar) !this unit is for exporting grid
       !if file exists, remove it
       export_file_var = TRIM(export_path) //  'wind_speed_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(dtWindSpeed,dtGrid) == 0)) THEN
            CALL Catch ('error', 'WindFlux',   &
                'dt wind speed is not multiple of dt of input grid')
        END IF
   ELSE
        !populate  metadata
        CALL ReadMetadata (fileunit, anemometersWS)

        !check spatial reference system
        IF ( .NOT. anemometersWS % mapping == mask % grid_mapping)  THEN
            CALL Catch ('info', 'WindFlux',   &
		              'converting coordinate of stations')
            !convert stations' coordinate
            point1 % system = anemometersWS % mapping
            point2 % system = mask % grid_mapping
            anemometersWS % mapping = mask % grid_mapping
            DO i = 1, anemometersWS % countObs
              point1 % easting = anemometersWS % obs (i) % xyz % easting
              point1 % northing = anemometersWS % obs (i) % xyz % northing
              point1 % elevation = anemometersWS % obs (i) % z
              CALL Convert (point1, point2)
              anemometersWS % obs (i) % xyz % easting = point2 % easting
              anemometersWS % obs (i) % xyz % northing = point2 % northing 
            END DO
        END IF
        
       
  END IF
   
RETURN
END SUBROUTINE WindFluxInit

!==============================================================================
!| Description:
!   Read wind speed data
SUBROUTINE WindFluxRead &
!
( 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 + - (dtWindSpeed - anemometersWS % timeIncrement)
   timeString = time_toread
   CALL Catch ('info', 'WindFlux',   &
		                 'read new wind speed data: ', &
                     argument = timeString)


   SELECT CASE (interpolationMethod)
     CASE (0) !input grid
       
          CALL ReadField (fileGrid,  time_toread, &
                        dtWindSpeed, dtGrid, &
                        'M', windSpeed, &
                         varName = windSpeed % var_name)
      
      
     CASE DEFAULT !requires interpolation
      !read new  wind speed data
      CALL ReadData (network = anemometersWS, fileunit = fileunit, &
                      time = time_toread, aggr_time = dtWindSpeed, &
                      aggr_type = 'ave', tresh = valid_prcn)
      
      IF (elevationdrift == 1) THEN !load wind direction
          CALL ReadData (network = anemometersWD, fileunit = fileunitWD, &
                      time = time_toread, aggr_time = dtWindSpeed, &
                      aggr_type = 'ave', tresh = valid_prcn)
      END IF
      
       !interpolate
        IF (interpolationMethod_assignment == 1) THEN !only one method is applied
            
            SELECT CASE (interpolationMethod) 
                CASE (1:3)
                
                CALL Interpolate (network = anemometersWS, &
                            method = interpolationMethod, &
                            near = neighbors, &
                            idw_power = idw_power, & 
                            anisotropy = krige_anisotropy, &
                            varmodel = krige_varmodel, &
                            lags = krige_lags, &
                            maxlag = krige_maxlag, &
                            grid = windSpeed, &
                            devst = grid_devst) 
                
                CASE (4)
                    
                    CALL WindMicromet (anemometersWS, anemometersWD, slope, &
                                       curvature, slope_az, micrometSlopeWT, &
                                       micrometCurvatureWT, windSpeed  )
                    
                CASE (5)
                    
                    CALL WindGonzalezLongatt (anemometersWS, anemometersWD, &
                                              dem, windSpeed)
                
            END SELECT
                
        ELSE
            !loop trough interpolation methods
            DO i = 1, 5
                IF (interpolationMethod_vector(i) > 0) THEN
                    SELECT CASE (interpolationMethod_vector(i))
                        CASE (1:3)
                            CALL Interpolate (network = anemometersWS, &
                                    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) 
                            
                        CASE (4)
                             CALL WindMicromet (anemometersWS, anemometersWD, slope, &
                                       curvature, slope_az, micrometSlopeWT, &
                                       micrometCurvatureWT, interpolatedMap (i) )
                        CASE (5)
                             CALL WindGonzalezLongatt (anemometersWS, anemometersWD, &
                                              dem, interpolatedMap (i) )
                    END SELECT
                    
                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
                        windSpeed % 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 SELECT

  !apply scale factor and offset, if defined
	IF (offset_value /= MISSING_DEF_REAL) THEN
	   CALL Offset (windSpeed, offset_value)
	END IF
	
	
	IF (scale_factor /= MISSING_DEF_REAL) THEN
	   CALL Scale (windSpeed, 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 = windSpeed % reference_time
        IF (needConversion) THEN
           CALL GridConvert (windSpeed, exportedGrid)
           CALL GridConvert (grid_devst, exportedGridVar)
        ELSE
           exportedGrid = windSpeed
           exportedGridVar = grid_devst
        END IF 

        SELECT CASE (export_format)
        CASE (1) !esri-ascii
              export_file = TRIM(export_path) //  TRIM(timeString) // &
                           '_windspeed.asc'
              CALL Catch ('info', 'WindFlux',   &
		                       '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) // &
                           '_windspeed_variance.asc'
                    CALL Catch ('info', 'WindFlux',   &
		                       '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) // &
                           '_windspeed'
              CALL Catch ('info', 'WindFlux',   &
		                       '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) // &
                               '_windspeed_variance'
                   CALL Catch ('info', 'WindFlux',   &
		                       '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', 'WindFlux',   &
                           'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGrid, export_file, 'wind_speed', 'append')
              
              IF (krige_var == 1) THEN !export kriging variance
                  CALL SetCurrentTime (time, exportedGridVar)
                  CALL Catch ('info', 'WindFlux',   &
		                       'exporting grid to file: ', &
                              argument = TRIM(export_file_var))
                 CALL ExportGrid (exportedGridVar, export_file_var, 'wind_speed_variance', 'append')
              END IF
        END SELECT
        timeNewExport = timeNewExport + export_dt
    END IF
  ENDIF

  !time forward
  timeNew = timeNew + dtWindSpeed
END IF

RETURN
END SUBROUTINE WindFluxRead


!==============================================================================
!| 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 = 'wind-speed') ) THEN
	       neighbors = IniReadInt ('nearest-points', ini, &
                                   section = 'wind-speed')
      ELSE !nearest-points missing
          CALL Catch ('error', 'WindFlux',   &
		                 '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 = 'wind-speed') ) THEN
	       neighbors = IniReadInt ('nearest-points', ini, &
                                   section = 'wind-speed')
           CALL Catch ('info', 'WindFlux', 'neighbors set to: ', argument = ToString (neighbors) )
      ELSE !nearest-points missing
          CALL Catch ('error', 'WindFlux',   &
		                 'nearest-points missing in meteo configuration file')
      END IF
      
      !kriging-variance
      IF (KeyIsPresent ('kriging-variance', ini, &
                         section = 'wind-speed') ) THEN
	       krige_var = IniReadInt ('kriging-variance', ini, &
                                   section = 'wind-speed')
           
           CALL Catch ('info', 'WindFlux',   &
		                 'kriging variance export set to ', argument = ToString(krige_var) )
      ELSE !kriging-variance missing
           krige_var = 0                                     
           CALL Catch ('info', 'WindFlux',   &
		                 'kriging variance export set to default (0)')
      END IF
      
      !kriging-anisotropy
      IF (KeyIsPresent ('kriging-anisotropy', ini, &
                         section = 'wind-speed') ) THEN
	       krige_anisotropy = IniReadInt ('kriging-anisotropy', ini, &
                                   section = 'wind-speed')
           CALL Catch ('info', 'WindFlux', 'kriging anisotropy set to: ', argument = ToString (krige_anisotropy) )
      ELSE !kriging-anisotropy missing
           krige_anisotropy = 0                                     
           CALL Catch ('info', 'WindFlux',   &
		                 'kriging anisotropy set to default (0)')
      END IF
      
      !kriging-variogram-model
      IF (KeyIsPresent ('kriging-variogram-model', ini, &
                         section = 'wind-speed') ) THEN
	       krige_varmodel = IniReadInt ('kriging-variogram-model', ini, &
                                   section = 'wind-speed')
           CALL Catch ('info', 'WindFlux', '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', 'WindFlux',   &
		                 'kriging variogram model set to default (2 exponential)')
      END IF
      
      !kriging-lags
      IF (KeyIsPresent ('kriging-lags', ini, &
                         section = 'wind-speed') ) THEN
	       krige_lags = IniReadInt ('kriging-lags', ini, &
                                   section = 'wind-speed')
           IF (krige_lags == 0) krige_lags = 15
           CALL Catch ('info', 'WindFlux', 'kriging lags set to: ', argument = ToString (krige_lags) )
      ELSE !kriging-variogram-model
           krige_lags = 15                                     
           CALL Catch ('info', 'WindFlux',   &
		                 'kriging lags set to default (15)')
      END IF
      
       !kriging-maxlag
      IF (KeyIsPresent ('kriging-maxlag', ini, &
                         section = 'wind-speed') ) THEN
	       krige_maxlag = IniReadInt ('kriging-maxlag', ini, &
                                   section = 'wind-speed')
           CALL Catch ('info', 'WindFlux', 'kriging max lag set to: ', argument = ToString (krige_maxlag) )
      ELSE !kriging-maxlag missing
           krige_maxlag = 0
          CALL Catch ('error', 'WindFlux',   &
		                 'kriging max lag set to default (0)')
      END IF
      
      
    CASE (4) !micromet
      CALL NewGrid (interpolatedMap (4), mask, 0.)
      
    CASE (5) !gonzalez
      CALL NewGrid (interpolatedMap (5), mask, 0.)
	
END SELECT


RETURN
END SUBROUTINE SetSpecificProperties


END MODULE WindFlux