WindFlux.f90 Source File
Manage wind speed data with arbitrary time cumulation
Source Code
!! 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