AirTemperatureDailyMean.f90 Source File
Manage mean daily air temperature
Source Code
!! 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