!----- LGPL -------------------------------------------------------------------- ! ! Copyright (C) Stichting Deltares, 2011-2013. ! ! This library is free software; you can redistribute it and/or ! modify it under the terms of the GNU Lesser General Public ! License as published by the Free Software Foundation version 2.1. ! ! This library is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! Lesser General Public License for more details. ! ! You should have received a copy of the GNU Lesser General Public ! License along with this library; if not, see . ! ! contact: delft3d.support@deltares.nl ! Stichting Deltares ! P.O. Box 177 ! 2600 MH Delft, The Netherlands ! ! All indications and logos of, and references to, "Delft3D" and "Deltares" ! are registered trademarks of Stichting Deltares, and remain the property of ! Stichting Deltares. All rights reserved. ! $Id$ ! $HeadURL$ !> This module instantiates FileReaders and their source Items, just as a GUI/framework instantiates kernels and their target Items. !! @author adri.mourits@deltares.nl !! @author stef.hummel@deltares.nl !! @author edwin.bos@deltares.nl module m_ec_provider use m_ec_typedefs use m_ec_parameters use m_ec_support use m_ec_filereader_read use m_ec_stringbuffer use m_ec_bcreader use time_module use m_ec_module use m_ec_message use m_ec_parameters use precision use string_module use netcdf use multi_file_io implicit none private public :: ecSetFileReaderProperties interface ecSetFileReaderProperties module procedure ecProviderInitializeFileReader end interface ecSetFileReaderProperties public :: ecAtLeastOnPointIsCorrection logical :: ecAtLeastOnPointIsCorrection = .false. contains ! ======================================================================= !> Initialize a new BCBlock item, which in turn constructs and initializes a filereader function ecProviderInitializeBCBlock(instancePtr, bcBlockId, yyyymmdd, tzone, fileReaderId, fileName, quantityName, plilabel, istat) result(success) logical :: success !< function status type(tEcInstance), pointer :: instancePtr !< intent(in) integer, intent(in) :: bcBlockId !< unique bcBlock id integer, intent(in) :: yyyymmdd !< kernel ref date real(hp), intent(in) :: tzone !< kernel time zone integer, intent(out) :: fileReaderId !< unique fileReader id character(*), intent(in) :: fileName !< relative path of data file character(*), intent(in) :: quantityName !< name of quantity, needed for structured input files (NetCDF and BC) character(*), intent(in) :: plilabel !< identify a (set of) pli-points integer, intent(out) :: istat !< Detailed result status. \see{m_ec_parameters}. type(tEcBCBlock), pointer :: bcBlockPtr !< BCBlock corresponding to bcBlockId type(tEcFileReader), pointer :: fileReaderPtr !< FileReader associated with the BC instance integer(kind=8) :: filehandle integer :: iostat, istat_ integer :: unit real(hp) :: ref_date success = .false. istat = EC_UNKNOWN_ERROR bcBlockPtr => ecSupportFindBCBlock(instancePtr, bcBlockId) if (.not.associated(bcBlockPtr)) return if (.not.ecSupportOpenExistingFileGnu(filehandle, fileName)) then istat = EC_DATA_NOTFOUND return end if if (.not.ecInitBC (filename, filehandle, quantityName, plilabel, bcBlockPtr, iostat)) return ! Every BC block (instance) needs an associated filereader referring to it fileReaderId = ecCreateFileReader(instancePtr) fileReaderPtr => ecSupportFindFileReader(instancePtr, fileReaderId) if (.not.associated(fileReaderPtr)) return if (bcBlockPtr%func==BC_FUNC_TSERIES) then if (.not.ecSupportTimestringToUnitAndRefdate(bcBlockPtr%timeunit, unit, ref_date)) then return endif ! Parsing the time string failed endif ! Timeseries function fileReaderPtr%bc => bcBlockPtr fileReaderPtr%ofType = provFile_bc if (.not.ecProviderInitializeTimeFrame(fileReaderPtr, yyyymmdd, fileReaderPtr%tframe%k_timezone)) return success = .false. select case(fileReaderPtr%bc%func) case (BC_FUNC_TSERIES) success = ecProviderCreateUniformItems(instancePtr, fileReaderPtr) case (BC_FUNC_HARMONIC, BC_FUNC_ASTRO) success = ecProviderCreateFourierItems(instancePtr, fileReaderPtr) case (BC_FUNC_HARMOCORR, BC_FUNC_ASTROCORR) !RL: Create filereader for corrections ? success = ecProviderAddCorrectionToFourierItems(instancePtr, fileReaderPtr) case (BC_FUNC_QHTABLE) success = ecProviderCreateQhtableItems(instancePtr, fileReaderPtr) case (BC_FUNC_TIM3D) success = ecProviderCreatet3DItems(instancePtr, fileReaderPtr) case default call setECMessage("ERROR: unknown function type.") end select end function ecProviderInitializeBCBlock ! ======================================================================= !> Initialize a new FileReader, by constructing the complete tree of source Items. !! On the opposite end of the EC-module is a kernel, which constructs the complete tree of target Items. recursive function ecProviderInitializeFileReader(instancePtr, fileReaderId, fileType, fileName, refdat, tzone, quantityName, forcingFile, bcBlockId) result(success) logical :: success !< function status type(tEcInstance), pointer :: instancePtr !< intent(in) integer, intent(in) :: fileReaderId !< unique FileReader id integer, intent(in) :: fileType !< type of data file, see provFile enumeration character(len=*), intent(in) :: fileName !< relative path of data file integer, intent(in) :: refdat !< Kernel's reference date, format: Gregorian yyyymmdd real(kind=hp), intent(in) :: tzone !< Kernel's timezone. character(len=*), optional, intent(in) :: quantityName !< name of quantity, needed for structured input files (NetCDF and BC) character(len=*), optional, intent(in) :: forcingFile !< name of quantity, needed for structured input files (NetCDF and BC) integer, optional, intent(in) :: bcBlockId !< if this filereader needs to be connected to a BC header block ! type(tEcFileReader), pointer :: fileReaderPtr !< FileReader corresponding to fileReaderId character(maxFileNameLen) :: fName !< relative path of data file, converted to the correct length integer :: i !< loop counter character(maxNameLen) :: l_quantityName !< explicit length version of quantityName integer :: iostat !< status returned from various file operations ! success = .false. fileReaderPtr => null() fName = '' l_quantityName = '' ! if (len_trim(fileName) > maxFileNameLen) then call setECMessage("ERROR: ec_provider::ecProviderInitializeFileReader: The filename string is too long.") return end if ! fName = fileName fileReaderPtr => ecSupportFindFileReader(instancePtr, fileReaderId) if (associated(fileReaderPtr)) then fileReaderPtr%ofType = fileType fileReaderPtr%fileName = fName fileReaderPtr%fileHandle = -1 ! The filereader itself has now an invalid filehandle if (fileReaderPtr%ofType == provFile_bc) then ! The argument bcBlockID refers to the ID of the bcBlock instance to be connected here, ! compulsory for this type of filereader stop 'HIER HEB IK ECHT NIKS MEER TE ZOEKEN' else success = ecSupportOpenExistingFile(fileReaderPtr%fileHandle, fileReaderPtr%fileName) success = .false. endif if (.not.ecProviderInitializeTimeFrame(fileReaderPtr, refdat, tzone)) return fileReaderPtr%nItems = 0 ! Create source Items and their contained types, based on file type and file header. if (present(quantityName)) then l_quantityName = quantityName if (.not.ecProviderCreateItems(instancePtr, fileReaderPtr, forcingFile, l_quantityName)) return else if (.not.ecProviderCreateItems(instancePtr, fileReaderPtr)) return end if end if success = .true. end function ecProviderInitializeFileReader ! ======================================================================= !> Create source Items and their contained types, based on file type and file header. function ecProviderCreateItems(instancePtr, fileReaderPtr, bctfilename, quantityname) result(success) logical :: success !< function status type(tEcInstance), pointer :: instancePtr !< intent(in) type(tEcFileReader), pointer :: fileReaderPtr !< intent(inout) character(maxNameLen), optional :: quantityname !< Names of the quantities read from file, needed for structured files (NetCDF), !< but also for bct-file character(maxNameLen), optional :: bctfilename !< file name of bct-file with data ! success = .false. select case(fileReaderPtr%ofType) case (provFile_undefined) call setECMessage("ERROR: ec_provider::ecProviderCreateItems: Unsupported file type.") case (provFile_uniform, provFile_unimagdir) success = ecProviderCreateUniformItems(instancePtr, fileReaderPtr) case (provFile_svwp) call setECMessage("ERROR: ec_provider::ecProviderCreateItems: Unsupported file type.") case (provFile_svwp_weight) call setECMessage("ERROR: ec_provider::ecProviderCreateItems: Unsupported file type.") case (provFile_arcinfo) success = ecProviderCreateArcinfoItems(instancePtr, fileReaderPtr) case (provFile_spiderweb) success = ecProviderCreateSpiderwebItems(instancePtr, fileReaderPtr) case (provFile_curvi) success = ecProviderCreateCurviItems(instancePtr, fileReaderPtr) case (provFile_curvi_weight) call setECMessage("ERROR: ec_provider::ecProviderCreateItems: Unsupported file type.") case (provFile_triangulation) call setECMessage("ERROR: ec_provider::ecProviderCreateItems: Unsupported file type.") case (provFile_triangulationmagdir) call setECMessage("ERROR: ec_provider::ecProviderCreateItems: Unsupported file type.") case (provFile_qhtable) success = ecProviderCreateQhtableItems(instancePtr, fileReaderPtr) case (provFile_poly_tim) if (present(bctfilename)) then if (.not.present(quantityname)) then call setECMessage("ERROR: ec_provider::ecProviderCreateItems: BC type, but no quantity name.") return endif success = ecProviderCreatePolyTimItemsBC(instancePtr, fileReaderPtr, bctfilename, quantityname) else success = ecProviderCreatePolyTimItems(instancePtr, fileReaderPtr) endif case (provFile_fourier) success = ecProviderCreateFourierItems(instancePtr, fileReaderPtr) case (provFile_grib) call setECMessage("ERROR: ec_provider::ecProviderCreateItems: Unsupported file type.") case (provFile_netcdf) if (present(quantityname)) then select case(quantityname) case ("rainfall") success = ecProviderCreateNetcdfItems(instancePtr, fileReaderPtr, quantityname) case default success = ecProviderCreateWaveNetcdfItems(instancePtr, fileReaderPtr, quantityname) end select else call setECMessage("ERROR: ec_provider::ecProviderCreateItems: NetCDF requires a quantity name.") end if case (provFile_t3D) success = ecProviderCreatet3DItems(instancePtr, fileReaderPtr) case default call setECMessage("ERROR: ec_provider::ecProviderCreateItems: Unknown file type.") end select end function ecProviderCreateItems ! ======================================================================= !> Create source Items and their contained types, from a qh-table file. function ecProviderCreateQhtableItems(instancePtr, fileReaderPtr) result(success) logical :: success !< function status type(tEcInstance), pointer :: instancePtr !< intent(in) type(tEcFileReader), pointer :: fileReaderPtr !< intent(inout) ! real(hp), dimension(:), allocatable :: discharges !< the table's discharge values real(hp), dimension(:), allocatable :: waterlevels !< the table's water level values integer :: nr_rows !< the number of rows in the table integer :: itemId !< helper variable integer :: quantityId !< helper variable integer :: elementSetId !< helper variable integer :: field0Id !< helper variable integer :: field1Id !< helper variable type(tEcItem), pointer :: item_discharge !< Item type(tEcItem), pointer :: item_waterlevel !< Item type(tEcItem), pointer :: item_slope !< Item type(tEcItem), pointer :: item_crossing !< Item integer :: i !< loop counter integer :: n1, n2 !< helper variables ! success = .true. item_discharge => null() item_waterlevel => null() item_slope => null() item_crossing => null() n1 = 12345 ! arbitrary large number ! ! Determine the number of rows and read the data. success = ecQhtableReadAll(fileReaderPtr, discharges, waterlevels, nr_rows) if (.not. success) return ! Create the item 'discharge'. quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'discharge'))) then success = .false. end if elementSetId = ecCreateElementSet(instancePtr) if (.not. (ecSetElementSetType(instancePtr, elementSetId, elmSetType_cartesian) .and. & ecSetElementSetNumberOfCoordinates(instancePtr, elementSetId, nr_rows))) then success = .false. end if field0Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field0Id, nr_rows))) then success = .false. end if field1Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field1Id, nr_rows))) then success = .false. end if itemId = ecCreateItem(instancePtr) if (.not. (ecSetItemRole(instancePtr, itemId, itemType_source) .and. & ecSetItemType(instancePtr, itemId, accessType_fileReader) .and. & ecSetItemQuantity(instancePtr, itemId, quantityId) .and. & ecSetItemElementSet(instancePtr, itemId, elementSetId) .and. & ecSetItemSourceT0Field(instancePtr, itemId, field0Id) .and. & ecSetItemSourceT1Field(instancePtr, itemId, field1Id))) then success = .false. else item_discharge => ecSupportFindItem(instancePtr, itemId) end if ! Create the item 'waterlevel'. quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'waterlevel'))) then success = .false. end if elementSetId = ecCreateElementSet(instancePtr) if (.not. (ecSetElementSetType(instancePtr, elementSetId, elmSetType_cartesian) .and. & ecSetElementSetNumberOfCoordinates(instancePtr, elementSetId, nr_rows))) then success = .false. end if field0Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field0Id, nr_rows))) then success = .false. end if field1Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field1Id, nr_rows))) then success = .false. end if itemId = ecCreateItem(instancePtr) if (.not. (ecSetItemRole(instancePtr, itemId, itemType_source) .and. & ecSetItemType(instancePtr, itemId, accessType_fileReader) .and. & ecSetItemQuantity(instancePtr, itemId, quantityId) .and. & ecSetItemElementSet(instancePtr, itemId, elementSetId) .and. & ecSetItemSourceT0Field(instancePtr, itemId, field0Id) .and. & ecSetItemSourceT1Field(instancePtr, itemId, field1Id))) then success = .false. else item_waterlevel => ecSupportFindItem(instancePtr, itemId) end if ! Create the item 'slope'. quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'slope'))) then success = .false. end if elementSetId = ecCreateElementSet(instancePtr) if (.not. (ecSetElementSetType(instancePtr, elementSetId, elmSetType_cartesian) .and. & ecSetElementSetNumberOfCoordinates(instancePtr, elementSetId, nr_rows))) then success = .false. end if field0Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field0Id, nr_rows))) then success = .false. end if field1Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field1Id, nr_rows))) then success = .false. end if itemId = ecCreateItem(instancePtr) if (.not. (ecSetItemRole(instancePtr, itemId, itemType_source) .and. & ecSetItemType(instancePtr, itemId, accessType_fileReader) .and. & ecSetItemQuantity(instancePtr, itemId, quantityId) .and. & ecSetItemElementSet(instancePtr, itemId, elementSetId) .and. & ecSetItemSourceT0Field(instancePtr, itemId, field0Id) .and. & ecSetItemSourceT1Field(instancePtr, itemId, field1Id))) then success = .false. else item_slope => ecSupportFindItem(instancePtr, itemId) end if ! Create the item 'crossing'. quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'crossing'))) then success = .false. end if elementSetId = ecCreateElementSet(instancePtr) if (.not. (ecSetElementSetType(instancePtr, elementSetId, elmSetType_cartesian) .and. & ecSetElementSetNumberOfCoordinates(instancePtr, elementSetId, nr_rows))) then success = .false. end if field0Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field0Id, nr_rows))) then success = .false. end if field1Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field1Id, nr_rows))) then success = .false. end if itemId = ecCreateItem(instancePtr) if (.not. (ecSetItemRole(instancePtr, itemId, itemType_source) .and. & ecSetItemType(instancePtr, itemId, accessType_fileReader) .and. & ecSetItemQuantity(instancePtr, itemId, quantityId) .and. & ecSetItemElementSet(instancePtr, itemId, elementSetId) .and. & ecSetItemSourceT0Field(instancePtr, itemId, field0Id) .and. & ecSetItemSourceT1Field(instancePtr, itemId, field1Id))) then success = .false. else item_crossing => ecSupportFindItem(instancePtr, itemId) end if ! ===== finish initialization of Fields ===== ! This FileReader has constant values, so we store them in T0. ! The values in T1 are never used, so to comply with the uniform EC-module mechanism, it is set to a large, increasable value. ! Rewind the file if (success) then item_slope%sourceT0FieldPtr%arr1dPtr = 0.0_hp item_crossing%sourceT0FieldPtr%arr1dPtr = 0.0_hp do i=1, nr_rows ! initialize Field T0 item_discharge%sourceT0FieldPtr%arr1dPtr(i) = discharges(i) item_waterlevel%sourceT0FieldPtr%arr1dPtr(i) = waterlevels(i) n2 = i-1 if (i > 1) then if (discharges(i) <= discharges(n2)) then call setECMessage("ERROR: ec_provider::ecProviderCreateQhtableItems: First column should be ordered increasingly.") success = .false. exit end if item_slope%sourceT0FieldPtr%arr1dPtr(n2) = (waterlevels(i)-waterlevels(n2))/(discharges(i)-discharges(n2)) item_crossing%sourceT0FieldPtr%arr1dPtr(n2) = waterlevels(n2) - item_slope%sourceT0FieldPtr%arr1dPtr(n2) * discharges(n2) end if end do item_discharge%sourceT0FieldPtr%timesteps = 0.0_hp item_waterlevel%sourceT0FieldPtr%timesteps = 0.0_hp item_slope%sourceT0FieldPtr%timesteps = 0.0_hp item_crossing%sourceT0FieldPtr%timesteps = 0.0_hp ! initialize Field T1 (arbitrary far future) item_discharge%sourceT1FieldPtr%timesteps = 10000.0_hp item_waterlevel%sourceT1FieldPtr%timesteps = 10000.0_hp item_slope%sourceT1FieldPtr%timesteps = 10000.0_hp item_crossing%sourceT1FieldPtr%timesteps = 10000.0_hp end if ! Add successfully created source Items to the FileReader if (success) success = ecAddItemToFileReader(instancePtr, fileReaderPtr%id, item_discharge%id) if (success) success = ecAddItemToFileReader(instancePtr, fileReaderPtr%id, item_waterlevel%id) if (success) success = ecAddItemToFileReader(instancePtr, fileReaderPtr%id, item_slope%id) if (success) success = ecAddItemToFileReader(instancePtr, fileReaderPtr%id, item_crossing%id) ! finalize if (allocated(discharges)) deallocate(discharges) if (allocated(waterlevels)) deallocate(waterlevels) end function ecProviderCreateQhtableItems ! ======================================================================= !> Create source Items and their contained types, based on Fourier file header. function ecProviderCreateFourierItems(instancePtr, fileReaderPtr) result(success) logical :: success !< function status type(tEcInstance), pointer :: instancePtr !< intent(in) type(tEcFileReader), pointer :: fileReaderPtr !< intent(inout) ! integer :: i !< loop counter integer :: quantityId !< helper variable integer :: elementSetId !< helper variable integer :: field0Id !< helper variable integer :: field1Id !< helper variable integer :: itemId !< helper variable ! integer :: nPeriods !< number of periods type(tEcItem), pointer :: itemPeriod !< Item containing the Fourier period type(tEcItem), pointer :: itemMagnitude !< Item containing the Fourier magnitude type(tEcItem), pointer :: itemPhase !< Item containing the Fourier phase real(hp), dimension(:), allocatable :: periods !< Fourier components transformed into periods real(hp), dimension(:), allocatable :: magnitudes !< seed values for the magnitudes of the Fourier components real(hp), dimension(:), allocatable :: phases !< seed values for the phases of the Fourier components integer :: yyyymmdd !< TimeFrame's reference date as Gregorian yyyymmdd ! success = .true. itemPeriod => null() itemMagnitude => null() itemPhase => null() ! ! Boundary condition components file containing: either period[minute] or component name, amplitude[m], phase[degree]. call jul2ymd(int(fileReaderPtr%tframe%k_refdate + 2400000.5_hp), yyyymmdd) if (ecFourierReadAll(fileReaderPtr, periods, magnitudes, phases, nPeriods, yyyymmdd)) then ! ===== all periods ===== quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'period') .and. & ecSetQuantityUnits(instancePtr, quantityId, 'minute'))) then success = .false. end if elementSetId = ecCreateElementSet(instancePtr) if (.not. (ecSetElementSetType(instancePtr, elementSetId, elmSetType_cartesian) .and. & ecSetElementSetNumberOfCoordinates(instancePtr, elementSetId, nPeriods))) then success = .false. end if field0Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field0Id, nPeriods))) then success = .false. end if field1Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field1Id, nPeriods))) then success = .false. end if itemId = ecCreateItem(instancePtr) if (.not. (ecSetItemRole(instancePtr, itemId, itemType_source) .and. & ecSetItemType(instancePtr, itemId, accessType_fileReader) .and. & ecSetItemQuantity(instancePtr, itemId, quantityId) .and. & ecSetItemElementSet(instancePtr, itemId, elementSetId) .and. & ecSetItemSourceT0Field(instancePtr, itemId, field0Id) .and. & ecSetItemSourceT1Field(instancePtr, itemId, field1Id))) then success = .false. else itemPeriod => ecSupportFindItem(instancePtr, itemId) end if ! ===== magnitude ===== quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'magnitude') .and. & ecSetQuantityUnits(instancePtr, quantityId, 'm'))) then success = .false. end if elementSetId = ecCreateElementSet(instancePtr) if (.not. (ecSetElementSetType(instancePtr, elementSetId, elmSetType_cartesian) .and. & ecSetElementSetNumberOfCoordinates(instancePtr, elementSetId, nPeriods))) then success = .false. end if field0Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field0Id, nPeriods))) then success = .false. end if field1Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field1Id, nPeriods))) then success = .false. end if itemId = ecCreateItem(instancePtr) if (.not. (ecSetItemRole(instancePtr, itemId, itemType_source) .and. & ecSetItemType(instancePtr, itemId, accessType_fileReader) .and. & ecSetItemQuantity(instancePtr, itemId, quantityId) .and. & ecSetItemElementSet(instancePtr, itemId, elementSetId) .and. & ecSetItemSourceT0Field(instancePtr, itemId, field0Id) .and. & ecSetItemSourceT1Field(instancePtr, itemId, field1Id))) then success = .false. else itemMagnitude => ecSupportFindItem(instancePtr, itemId) end if ! ===== phase ===== quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'phase') .and. & ecSetQuantityUnits(instancePtr, quantityId, 'degree'))) then success = .false. end if elementSetId = ecCreateElementSet(instancePtr) if (.not. (ecSetElementSetType(instancePtr, elementSetId, elmSetType_cartesian) .and. & ecSetElementSetNumberOfCoordinates(instancePtr, elementSetId, nPeriods))) then success = .false. end if field0Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field0Id, nPeriods))) then success = .false. end if field1Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field1Id, nPeriods))) then success = .false. end if itemId = ecCreateItem(instancePtr) if (.not. (ecSetItemRole(instancePtr, itemId, itemType_source) .and. & ecSetItemType(instancePtr, itemId, accessType_fileReader) .and. & ecSetItemQuantity(instancePtr, itemId, quantityId) .and. & ecSetItemElementSet(instancePtr, itemId, elementSetId) .and. & ecSetItemSourceT0Field(instancePtr, itemId, field0Id) .and. & ecSetItemSourceT1Field(instancePtr, itemId, field1Id))) then success = .false. else itemPhase => ecSupportFindItem(instancePtr, itemId) end if ! ===== finish initialization of Fields ===== ! This FileReader calculates it's values, rather then reading them in, so we store the seed values in T0. ! The values in T1 are never used, so to comply with the uniform EC-module mechanism, it is set to a large, increasable value. ! The actual calculation is performed by the Converter. if (success) then do i=1, nPeriods ! initialize Field T0 itemPeriod%sourceT0FieldPtr%arr1dPtr(i) = periods(i) itemMagnitude%sourceT0FieldPtr%arr1dPtr(i) = magnitudes(i) itemPhase%sourceT0FieldPtr%arr1dPtr(i) = phases(i) end do itemPeriod%sourceT0FieldPtr%timesteps = 0.0_hp itemMagnitude%sourceT0FieldPtr%timesteps = 0.0_hp itemPhase%sourceT0FieldPtr%timesteps = 0.0_hp ! initialize Field T1 (arbitrary far future) itemPeriod%sourceT1FieldPtr%timesteps = 10000.0_hp itemMagnitude%sourceT1FieldPtr%timesteps = 10000.0_hp itemPhase%sourceT1FieldPtr%timesteps = 10000.0_hp end if ! Add successfully created source Items to the FileReader if (success) success = ecAddItemToFileReader(instancePtr, fileReaderPtr%id, itemPeriod%id) if (success) success = ecAddItemToFileReader(instancePtr, fileReaderPtr%id, itemMagnitude%id) if (success) success = ecAddItemToFileReader(instancePtr, fileReaderPtr%id, itemPhase%id) else success = .false. end if end function ecProviderCreateFourierItems ! ======================================================================= !> Add corr. to astro/harmonic components function ecProviderAddCorrectionToFourierItems(instancePtr, corrFileReaderPtr) result(success) logical :: success !< function status type(tEcInstance), pointer :: instancePtr !< intent(in) type(tEcFileReader), pointer :: corrFileReaderPtr !< intent(inout) ! integer :: i !< loop counter integer :: quantityId !< helper variable integer :: elementSetId !< helper variable integer :: field0Id !< helper variable integer :: field1Id !< helper variable integer :: itemId !< helper variable ! integer :: nPeriods !< number of periods real(hp), dimension(:), allocatable :: periods !< Fourier components transformed into periods real(hp), dimension(:), allocatable :: magnitudes !< seed values for the magnitudes of the Fourier components real(hp), dimension(:), allocatable :: phases !< seed values for the phases of the Fourier components integer :: yyyymmdd !< TimeFrame's reference date as Gregorian yyyymmdd type(tEcFileReader) , pointer :: cmpFileReaderPtr !< related file reader (with components) type (tEcItem), pointer :: cmp_period, cmp_amplitude, cmp_phase real(hp), pointer :: cmpperiods(:), cmpamplitudes(:), cmpphases(:) integer :: icmp, ncmp, icor, ncor, iitem logical :: cmpfound ! success = .true. cmp_period => null() cmp_amplitude => null() cmp_phase => null() if (corrFileReaderPtr%bc%func==BC_FUNC_HARMOCORR) then cmpFileReaderPtr => ecProviderFindRelatedBCBlock(instancePtr, corrFileReaderPtr, corrFileReaderPtr % bc % qname, corrFileReaderPtr % bc % bcname, BC_FUNC_HARMONIC) elseif (corrFileReaderPtr%bc%func==BC_FUNC_ASTROCORR) then cmpFileReaderPtr => ecProviderFindRelatedBCBlock(instancePtr, corrFileReaderPtr, corrFileReaderPtr % bc % qname, corrFileReaderPtr % bc % bcname, BC_FUNC_ASTRO) endif if (.not. associated(cmpFileReaderPtr)) then ! TODO: message: no related component success = .false. return endif call jul2ymd(int(cmpFileReaderPtr%tframe%k_refdate + 2400000.5_hp), yyyymmdd) if (ecFourierReadAll(corrFileReaderPtr, periods, magnitudes, phases, nPeriods, yyyymmdd)) then do iitem = 1, cmpFileReaderPtr%nItems if (cmpFileReaderPtr%items(iitem)%ptr%role == itemType_source) then ! source items select case (cmpFileReaderPtr%items(iitem)%ptr%quantityptr%name) case('period') cmp_period => cmpFileReaderPtr%items(iitem)%ptr case('magnitude') cmp_amplitude => cmpFileReaderPtr%items(iitem)%ptr case('phase') cmp_phase => cmpFileReaderPtr%items(iitem)%ptr end select endif enddo ncmp = size(cmp_period%sourceT0FieldPtr%arr1d) ! number of components ncor = size(periods) ! number of corrected components cmpperiods => cmp_period%sourceT0FieldPtr%arr1d cmpamplitudes => cmp_amplitude%sourceT0FieldPtr%arr1d cmpphases => cmp_phase%sourceT0FieldPtr%arr1d do icor = 1, ncor cmpfound = .false. do icmp = 1, ncmp if (comparereal(periods(icor),cmpperiods(icmp))==0) then cmpamplitudes(icmp)=cmpamplitudes(icmp)*magnitudes(icor) cmpphases(icmp)=cmpphases(icmp)+phases(icor) cmpfound = .true. cycle endif enddo ! components if(.not.cmpfound) then ! TODO: think of a meaningfull error message if correcting a non-existing component success = .false. endif enddo ! corrections ! TODO: ! find 'periods' source item in cmpFileReaderPtr ! find 'magnitudes' source item in cmpFileReaderPtr ! find 'phases' source item in cmpFileReaderPtr ! loop over periods in correction ! periodFound = false ! for all periods in 'periods' source item ! if (period == period in source item) ! periodFound = true ! add phase, multiply magnitude ! cycle ! if (not periodFound) ERROR else ! TODO: message success = .false. end if end function ecProviderAddCorrectionToFourierItems ! ======================================================================= !> Find the file reader for the bc block that contains the component definition for the comp.correction block function ecProviderFindRelatedBCBlock(instancePtr, corrFileReaderPtr, qname, bcname, func) result(cmpFileReaderPtr) type(tEcFileReader), pointer :: cmpFileReaderPtr !< resulting file reader type(tEcInstance), pointer :: instancePtr !< intent(in) type(tEcFileReader), pointer :: corrFileReaderPtr !< intent(inout) character(len=*) :: qname !< quantity name character(len=*) :: bcname !< point on poly name integer, intent(in) :: func !< function type integer :: iFileReader type (tEcBCBlock), pointer :: BCBlockptr ! cmpFileReaderPtr => null() do iFileReader = 1, instancePtr%nFileReaders BCBlockptr => instancePtr%EcFileReadersPtr(iFileReader)%ptr%bc if (associated(BCBlockptr)) then if (trim(BCBlockptr%bcname)==trim(bcname) .and. (trim(BCBlockptr%qname)==trim(qname)) & .and. BCBlockptr%func == func) then cmpFileReaderPtr => instancePtr%EcFileReadersPtr(iFileReader)%ptr exit endif else endif enddo end function ecProviderFindRelatedBCBlock ! ======================================================================= !> Create source Item and its contained types, for timeseries source data, !! which may come from a time series file ('uniform'), or a BC-file block. !! Uniform file records (rows) are treated as vector quantities, reading all physical quantities (columns) into a single Item. function ecProviderCreateUniformItems(instancePtr, fileReaderPtr) result(success) logical :: success !< function status type(tEcInstance), pointer :: instancePtr !< intent(in) type(tEcFileReader), pointer :: fileReaderPtr !< intent(inout) ! character(len=132) :: rec !< first data line in file integer :: n_quantities !< number of quantities in the file integer :: quantityId !< helper variable integer :: elementSetId !< helper variable integer :: field0Id !< helper variable integer :: field1Id !< helper variable integer :: itemId !< helper variable integer :: indx !< real(hp) :: time_steps !< number of time steps for next data block type(tEcItem), pointer :: item !< Item containing all components character(len=1) :: excl !< ! success = .true. item => null() ! select case (fileReaderPtr%ofType) case (provFile_uniform, provFile_unimagdir) rec = ecUniReadFirstLine(fileReaderPtr) excl = '!' indx = index(rec, excl) if (indx /= 0) then n_quantities = count_words(rec(1:indx-1)) - 1 else n_quantities = count_words(rec) - 1 end if case (provFile_bc) n_quantities = count(fileReaderPtr%bc%quantity%jacolumn) end select ! time [minute], quantities quantityId = ecCreateQuantity(instancePtr) select case (fileReaderPtr%ofType) case (provFile_uniform, provFile_unimagdir) if (.not. ecSetQuantityName(instancePtr, quantityId, 'uniform_item')) then success = .false. end if case (provFile_bc) if (.not. ecSetQuantityName(instancePtr, quantityId, 'uniform_item')) then ! trim(fileReaderPtr%bc%qname))) then success = .false. end if end select ! N_quantities number of scalar quantities. elementSetId = ecCreateElementSet(instancePtr) if (.not. (ecSetElementSetType(instancePtr, elementSetId, elmSetType_scalar) .and. & ecSetElementSetNumberOfCoordinates(instancePtr, elementSetId, n_quantities))) then success = .false. end if ! N_quantities scalars in a Field array. field0Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field0Id, n_quantities))) then success = .false. end if field1Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field1Id, n_quantities))) then success = .false. end if itemId = ecCreateItem(instancePtr) if (.not. (ecSetItemRole(instancePtr, itemId, itemType_source) .and. & ecSetItemType(instancePtr, itemId, accessType_fileReader) .and. & ecSetItemQuantity(instancePtr, itemId, quantityId) .and. & ecSetItemElementSet(instancePtr, itemId, elementSetId) .and. & ecSetItemSourceT0Field(instancePtr, itemId, field0Id) .and. & ecSetItemSourceT1Field(instancePtr, itemId, field1Id))) then success = .false. else item => ecSupportFindItem(instancePtr, itemId) end if !select case (fileReaderPtr%ofType) ! case (provFile_uniform, provFile_unimagdir) ! if (.not. ecSetQuantityName(instancePtr, quantityId, 'uniform_item')) then ! success = .false. ! end if ! case (provFile_bc) ! if (.not. ecSetQuantityName(instancePtr, quantityId, trim(fileReaderPtr%bc%qname))) then ! success = .false. ! end if !end select ! ===== finish initialization of Fields ===== ! Read the first two records into tEcItem%sourceT0FieldPtr and tEcItem%sourceT1FieldPtr. select case (fileReaderPtr%ofType) case (provFile_uniform, provFile_unimagdir) rewind(unit=fileReaderPtr%fileHandle) if (success) then success = ecUniReadBlock(fileReaderPtr, item%sourceT0FieldPtr%timesteps, item%sourceT0FieldPtr%arr1dPtr) endif if (success) then success = ecUniReadBlock(fileReaderPtr, item%sourceT1FieldPtr%timesteps, item%sourceT1FieldPtr%arr1dPtr) endif case (provFile_bc) if (success) then success = ecBCReadBlock(fileReaderPtr, item%sourceT0FieldPtr%timesteps, item%sourceT0FieldPtr%arr1dPtr) endif if (success) then success = ecBCReadBlock(fileReaderPtr, item%sourceT1FieldPtr%timesteps, item%sourceT1FieldPtr%arr1dPtr) endif end select ! Add successfully created source Item to the FileReader if (success) success = ecAddItemToFileReader(instancePtr, fileReaderPtr%id, item%id) end function ecProviderCreateUniformItems ! ======================================================================= !> Create source Items and their contained types, based on Arcinfo file header. function ecProviderCreateArcinfoItems(instancePtr, fileReaderPtr) result(success) logical :: success !< function status type(tEcInstance), pointer :: instancePtr !< intent(in) type(tEcFileReader), pointer :: fileReaderPtr !< intent(inout) ! integer :: n_cols !< number of columns integer :: n_rows !< number of rows real(hp) :: x0 !< seed x-coordinate real(hp) :: y0 !< seed y-coordinate real(hp) :: dxa !< step size in x real(hp) :: dya !< step size in y real(hp) :: dmiss !< missing data value integer :: quantityId !< helper variable integer :: elementSetId !< helper variable integer :: field0Id !< helper variable integer :: field1Id !< helper variable integer :: itemId !< helper variable type(tEcItem), pointer :: item character(len=20) :: name !< integer :: i !< loop counter ! success = .true. item => null() ! Read relevant information from the first arcinfo file's header. if (readarcinfoheader(fileReaderPtr%fileHandle, n_cols, n_rows, x0, y0, dxa, dya, dmiss)) then ! One common ElementSet. elementSetId = ecCreateElementSet(instancePtr) if (.not. (ecSetElementSetType(instancePtr, elementSetId, elmSetType_cartesian_equidistant) .and. & ecSetElementSetX0Dx(instancePtr, elementSetId, x0, dxa) .and. & ecSetElementSetY0Dy(instancePtr, elementSetId, y0, dya) .and. & ecSetElementSetRowsCols(instancePtr, elementSetId, n_rows, n_cols))) then success = .false. end if ! if (index(fileReaderPtr%fileName, '.amu') /= 0) then ! ===== quantity: wind component u (usually == x) ===== quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'wind_u') .and. & ecSetQuantityUnits(instancePtr, quantityId, trim(ecSpiderwebAndCurviFindInFile(fileReaderPtr%fileHandle, 'unit1'))))) then success = .false. end if else if (index(fileReaderPtr%fileName, '.amv') /= 0) then ! ===== quantity: wind component v (usually == y) ===== quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'wind_v') .and. & ecSetQuantityUnits(instancePtr, quantityId, trim(ecSpiderwebAndCurviFindInFile(fileReaderPtr%fileHandle, 'unit1'))))) then success = .false. end if else if (index(fileReaderPtr%fileName, '.amp') /= 0) then ! ===== quantity: wind component p ===== quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'wind_p') .and. & ecSetQuantityUnits(instancePtr, quantityId, trim(ecSpiderwebAndCurviFindInFile(fileReaderPtr%fileHandle, 'unit1'))))) then success = .false. end if end if field0Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field0Id, n_cols*n_rows) .and. & ecSetFieldMissingValue(instancePtr, field0Id, dmiss))) then success = .false. end if field1Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field1Id, n_cols*n_rows) .and. & ecSetFieldMissingValue(instancePtr, field1Id, dmiss))) then success = .false. end if itemId = ecCreateItem(instancePtr) if (.not. ( ecSetItemRole(instancePtr, itemId, itemType_source) .and. & ecSetItemType(instancePtr, itemId, accessType_fileReader) .and. & ecSetItemQuantity(instancePtr, itemId, quantityId) .and. & ecSetItemElementSet(instancePtr, itemId, elementSetId) .and. & ecSetItemSourceT0Field(instancePtr, itemId, field0Id) .and. & ecSetItemSourceT1Field(instancePtr, itemId, field1Id))) then success = .false. else item => ecSupportFindItem(instancePtr, itemId) end if else success = .false. ! failed reading header end if ! ===== finish initialization of Fields ===== ! Read the first two records into tEcItem%sourceT0FieldPtr and tEcItem%sourceT1FieldPtr. if (success) then rewind(unit=fileReaderPtr%fileHandle) success = ecSpiderAndCurviAndArcinfoReadToBody(fileReaderPtr%fileHandle) end if if (success) then success = ecArcinfoAndCurviReadBlock(fileReaderPtr, fileReaderPtr%fileHandle, item, 0, n_cols, n_rows) end if if (success) then success = ecArcinfoAndCurviReadBlock(fileReaderPtr, fileReaderPtr%fileHandle, item, 1, n_cols, n_rows) end if ! Add successfully created source Item to the FileReader if (success) success = ecAddItemToFileReader(instancePtr, fileReaderPtr%id, item%id) end function ecProviderCreateArcinfoItems ! ======================================================================= !> Construct an ElementSet from a curvilinear grid file. !! meteo1: reaarc_curv_tim function ecProviderCreateCurviElementSet(instancePtr, fileReaderPtr, elementSetId, n_cols, n_rows) result(success) logical :: success !< function status type(tEcInstance), pointer :: instancePtr !< intent(in) type(tEcFileReader), pointer :: fileReaderPtr !< intent(inout) integer, intent(out) :: elementSetId !< if of new ElementSet integer, intent(out) :: n_cols, n_rows ! character(len=132) :: rec !< a read line character(len=maxFileNameLen) :: grid_file !< file name of curvilinear grid integer :: minp !< IO unit number integer :: mx, nx !< n_clos, n_rows real(hp), dimension(:), allocatable :: x, y !< coordinate arrays integer :: i, j !< loop counters character(len=10) :: dummy !< helper variable for ignored data integer :: istat !< status of operation integer :: elmSetType !< spherical (elmSetType_spheric) or not (elmSetType_Cartesian) ! success = .false. ! Find and open the curvilinear grid file. rec = ecSpiderwebAndCurviFindInFile(fileReaderPtr%fileHandle, 'grid_file') read(rec, *) grid_file success = ecSupportOpenExistingFile(minp, grid_file) if (.not. success) return ! Read the file header. elmSetType = elmSetType_Cartesian 20 read(minp,'(a)') rec if (index (rec,'=') == 0) then ! Backwards compatible: first line could contain spherical keyword if (index(rec, 'Spherical') >= 1 .or. & index(rec, 'SPHERICAL') >= 1 ) then ! grid has spherical coordinates. elmSetType = elmSetType_spheric endif goto 20 end if read(minp,*) mx, nx read(minp,*) ! skips a line ! Read the file body. allocate(x(mx*nx)) allocate(y(mx*nx)) ! Read the data row-by-row into a 1D array. do j=0, nx-1 read(minp, *, iostat = istat) dummy, dummy, (x(j*mx+i), i=1, mx) end do do j=0, nx-1 read(minp, *, iostat = istat) dummy, dummy, (y(j*mx+i), i=1, mx) end do close(minp) if (istat /= 0) then call setECMessage("ERROR: ec_provider::ecProviderCreateCurviElementSet: Unable to read data block from curvilinear grid file.") success = .false. return end if ! elementSetId = ecCreateElementSet(instancePtr) if (.not. (ecSetElementSetType(instancePtr, elementSetId, elmSetType) .and. & ecSetElementSetXArray(instancePtr, elementSetId, x) .and. & ecSetElementSetYArray(instancePtr, elementSetId, y) .and. & ecSetElementSetRowsCols(instancePtr, elementSetId, nx, mx))) then success = .false. end if deallocate(x) deallocate(y) n_cols = mx n_rows = nx end function ecProviderCreateCurviElementSet ! ======================================================================= !> Create source Items with their contained types, based on curvilinear file header. !! metteo1: reaarc_curv_tim function ecProviderCreateCurviItems(instancePtr, fileReaderPtr) result(success) logical :: success !< function status type(tEcInstance), pointer :: instancePtr !< intent(in) type(tEcFileReader), pointer :: fileReaderPtr !< intent(inout) ! integer :: quantityId !< if of new Quantity integer :: elementSetId !< if of new ElementSet integer :: field0Id !< if of new Field integer :: field1Id !< if of new Field integer :: itemId !< if of new Item type(tEcItem), pointer :: wind_p !< Item containing wind_p type(tEcItem), pointer :: wind_x !< Item containing wind_x type(tEcItem), pointer :: wind_y !< Item containing wind_y character(len=132) :: rec !< a read line real(hp) :: missingValue !< helper variable integer :: n_cols, n_rows integer :: check logical :: current_block !< indicates whether to read the next quantity from the current block or from the next block ! success = .false. wind_p => null() wind_x => null() wind_y => null() ! rec = ecSpiderwebAndCurviFindInFile(fileReaderPtr%fileHandle, 'NODATA_value') if (len_trim(rec) /= 0) then read(rec, *) missingValue else missingValue = -9999.0_hp end if ! rec = ecSpiderwebAndCurviFindInFile(fileReaderPtr%fileHandle, 'n_quantity') if (len_trim(rec) == 0) then call setECMessage("ERROR: ec_provider::ecProviderCreateCurviItems: Invalid file header specified.") return end if read(rec, *) check if (check /= 3) then call setECMessage("ERROR: ec_provider::ecProviderCreateCurviItems: Invalid number of quantities specified in header.") return end if ! success = ecProviderCreateCurviElementSet(instancePtr, fileReaderPtr, elementSetId, n_cols, n_rows) if (.not. success) return ! ! ===== quantity1: wind_p ===== quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'wind_p') .and. & ecSetQuantityUnits(instancePtr, quantityId, ' '))) then success = .false. end if field0Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field0Id, n_cols*n_rows) .and. & ecSetFieldMissingValue(instancePtr, field0Id, missingValue))) then success = .false. end if field1Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field1Id, n_cols*n_rows) .and. & ecSetFieldMissingValue(instancePtr, field1Id, missingValue))) then success = .false. end if itemId = ecCreateItem(instancePtr) if (.not. ( ecSetItemRole(instancePtr, itemId, itemType_source) .and. & ecSetItemType(instancePtr, itemId, accessType_fileReader) .and. & ecSetItemQuantity(instancePtr, itemId, quantityId) .and. & ecSetItemElementSet(instancePtr, itemId, elementSetId) .and. & ecSetItemSourceT0Field(instancePtr, itemId, field0Id) .and. & ecSetItemSourceT1Field(instancePtr, itemId, field1Id))) then success = .false. else wind_p => ecSupportFindItem(instancePtr, itemId) end if ! ===== quantity2: wind_x ===== quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'wind_x') .and. & ecSetQuantityUnits(instancePtr, quantityId, ' '))) then success = .false. end if field0Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field0Id, n_cols*n_rows) .and. & ecSetFieldMissingValue(instancePtr, field0Id, missingValue))) then success = .false. end if field1Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field1Id, n_cols*n_rows) .and. & ecSetFieldMissingValue(instancePtr, field1Id, missingValue))) then success = .false. end if itemId = ecCreateItem(instancePtr) if (.not. ( ecSetItemRole(instancePtr, itemId, itemType_source) .and. & ecSetItemType(instancePtr, itemId, accessType_fileReader) .and. & ecSetItemQuantity(instancePtr, itemId, quantityId) .and. & ecSetItemElementSet(instancePtr, itemId, elementSetId) .and. & ecSetItemSourceT0Field(instancePtr, itemId, field0Id) .and. & ecSetItemSourceT1Field(instancePtr, itemId, field1Id))) then success = .false. else wind_x => ecSupportFindItem(instancePtr, itemId) end if ! ===== quantity3: wind_y ===== quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'wind_y') .and. & ecSetQuantityUnits(instancePtr, quantityId, ' '))) then success = .false. end if field0Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field0Id, n_cols*n_rows) .and. & ecSetFieldMissingValue(instancePtr, field0Id, missingValue))) then success = .false. end if field1Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field1Id, n_cols*n_rows) .and. & ecSetFieldMissingValue(instancePtr, field1Id, missingValue))) then success = .false. end if itemId = ecCreateItem(instancePtr) if (.not. ( ecSetItemRole(instancePtr, itemId, itemType_source) .and. & ecSetItemType(instancePtr, itemId, accessType_fileReader) .and. & ecSetItemQuantity(instancePtr, itemId, quantityId) .and. & ecSetItemElementSet(instancePtr, itemId, elementSetId) .and. & ecSetItemSourceT0Field(instancePtr, itemId, field0Id) .and. & ecSetItemSourceT1Field(instancePtr, itemId, field1Id))) then success = .false. else wind_y => ecSupportFindItem(instancePtr, itemId) end if ! ===== finish initialization of Fields ===== ! Read the first two records into tEcItem%sourceT0FieldPtr and tEcItem%sourceT1FieldPtr. if (success) then rewind(unit=fileReaderPtr%fileHandle) success = ecSpiderAndCurviAndArcinfoReadToBody(fileReaderPtr%fileHandle) current_block = .true. end if if (success) success = ecArcinfoAndCurviReadBlock(fileReaderPtr, fileReaderPtr%fileHandle, wind_p, 0, n_cols, n_rows) if (success) success = ecArcinfoAndCurviReadBlock(fileReaderPtr, fileReaderPtr%fileHandle, wind_x, 0, n_cols, n_rows, current_block) if (success) success = ecArcinfoAndCurviReadBlock(fileReaderPtr, fileReaderPtr%fileHandle, wind_y, 0, n_cols, n_rows, current_block) if (success) success = ecArcinfoAndCurviReadBlock(fileReaderPtr, fileReaderPtr%fileHandle, wind_p, 1, n_cols, n_rows) if (success) success = ecArcinfoAndCurviReadBlock(fileReaderPtr, fileReaderPtr%fileHandle, wind_x, 1, n_cols, n_rows, current_block) if (success) success = ecArcinfoAndCurviReadBlock(fileReaderPtr, fileReaderPtr%fileHandle, wind_y, 1, n_cols, n_rows, current_block) ! Add successfully created source Items to the FileReader if (success) success = ecAddItemToFileReader(instancePtr, fileReaderPtr%id, wind_p%id) if (success) success = ecAddItemToFileReader(instancePtr, fileReaderPtr%id, wind_x%id) if (success) success = ecAddItemToFileReader(instancePtr, fileReaderPtr%id, wind_y%id) end function ecProviderCreateCurviItems ! ======================================================================= !> Create source Items with their contained types, based on t3D file header. function ecProviderCreatet3DItems(instancePtr, fileReaderPtr) result(success) logical :: success !< function status type(tEcInstance), pointer :: instancePtr !< intent(in) type(tEcFileReader), pointer :: fileReaderPtr !< intent(inout) ! integer :: quantityId !< id of new Quantity integer :: elementSetId !< id of new ElementSet integer :: field0Id !< id of new Field integer :: field1Id !< id of new Field integer :: itemId !< id of new Item integer :: layer_type !< type of layer type(tEcItem), pointer :: valueptr !< Item containing z/sigma-dependent values type(tEcBCBlock), pointer :: bcptr character(len=132) :: rec !< a read line integer, parameter :: MAXSTRLEN=128 !< integer, parameter :: MAXLAY=256 !< integer :: numlay, i !< real(hp), dimension(:), allocatable :: zws !< z-values of vertical velocities real(hp), dimension(MAXLAY) :: a !< integer :: ilay, jlay integer :: col_tmp real(hp) :: zws_tmp ! success = .false. valueptr => null() select case (fileReaderPtr%ofType) case (provFile_t3D) a = ec_undef_hp ! ! check if these are z-layers or sigma-layers rec = ecSpiderwebAndCurviFindInFile(fileReaderPtr%fileHandle, 'LAYER_TYPE') if (len_trim(rec) == 0) then call setECMessage("ERROR: ec_provider::ecProviderCreatet3DItems: Unable to find LAYER_TYPE in file header.") return end if ! call str_lower(rec) if ( index(rec,'sigma') /= 0 ) then layer_type = 0 else if ( index(rec,'z') /= 0 ) then layer_type = 1 else call setECMessage("ERROR: ec_provider::ecProviderCreatet3DItems: Invalid LAYER_TYPE specified in header.") return end if ! ! Acquire the number of layers by the number of readable columns in LAYERS field of header rec = ect3DFindInFile(fileReaderPtr%fileHandle, 'LAYERS') if (len_trim(rec) == 0) then call setECMessage("ERROR: ec_provider::ecProviderCreatet3DItems: Unable to find LAYER in file header.") return end if ! read(rec,*, end = 666) (a(i),i=1,MAXLAY) 666 continue ! numlay = 0 do i=1,MAXLAY if ( a(i) /= ec_undef_hp ) then numlay = numlay+1 end if end do ! ! allocate vertical mesh array if ( numlay > 0 ) then allocate(zws(numlay)) else call setECMessage("ERROR: ec_provider::ecProviderCreatet3DItems: no layers found in header") return end if ! zws = (/ (a(i), i=1,numlay) /) case (provFile_bc) ! Check if these are z-layers or sigma-layers if (.not.associated(fileReaderPtr%bc)) then call setECMessage("ERROR: ec_provider::ecProviderCreatet3DItems: bc-filetype, but no bc instance associated to filereader") return endif bcptr => fileReaderPtr%bc select case (bcptr%vptyp) case (BC_VPTYP_PERCBED) layer_type = 0 case (BC_VPTYP_ZBED) layer_type = 1 case default ! Invalid vertical position specification type call setECMessage("ERROR: ec_provider::ecProviderCreatet3DItems: Invalid vertical position specification type.") return end select ! Defined types (yet to be implemented): ! integer, parameter :: BC_VPTYP_PERCBED = 1 !< precentage from bed ! integer, parameter :: BC_VPTYP_ZDATUM = 2 !< z above datum ! integer, parameter :: BC_VPTYP_BEDSURF = 3 !< bedsurface ! integer, parameter :: BC_VPTYP_PERCSURF = 4 !< percentage from surface ! integer, parameter :: BC_VPTYP_ZBED = 5 !< z from bed ! integer, parameter :: BC_VPTYP_ZSURF = 6 !< z from surface ! Acquire the number of layers, ! Corresponds with the number of matching quantity blocks in a bc-header numlay = count(bcptr%quantity%jacolumn(1:bcptr%numcols)) ! Prepare the z-coordinates for the levels in zws(dimension=numlay) ! Corresponds with the SELECTED vertical positions in vp(:) ! allocate vertical mesh array if ( numlay > 0 ) then allocate(zws(numlay)) else call setECMessage("ERROR: ec_provider::ecProviderCreatet3DItems: no layers found in header") return end if zws = (/ (bcptr%vp(i), i=1,numlay) /) ! sort (zws,bcptr%quantity%vect2col) to later to be able to sort vector by vertical coordinate do jlay=numlay,1,-1 zws_tmp = zws(jlay) col_tmp = bcptr%quantity%vect2col(jlay) do ilay=jlay-1,1,-1 if (zws_tmp < zws(ilay)) then zws(ilay+1) = zws(ilay) bcptr%quantity%vect2col(ilay+1) = bcptr%quantity%vect2col(ilay) else zws(ilay+1) = zws_tmp bcptr%quantity%vect2col(ilay+1) = col_tmp exit endif enddo if (ilay==0) then zws(1) = zws_tmp bcptr%quantity%vect2col(1) = col_tmp endif enddo case default call setECMessage("ERROR: ec_provider::ecProviderCreatet3DItems: Unknown file type.") return end select ! ! ===== single quantity: quant ===== (Further code in this sub independent of file type) quantityId = ecCreateQuantity(instancePtr) if (.not.ecSetQuantityName(instancePtr, quantityId, 'quant')) return if (.not.ecSetQuantityUnits(instancePtr, quantityId, ' ')) return elementSetId = ecCreateElementSet(instancePtr) if (.not.ecSetElementSetType(instancePtr, elementSetId, elmSetType_cartesian)) return if (.not.ecSetElementSetXArray(instancePtr, elementSetId, zws)) return if (.not.ecSetElementSetNumberOfCoordinates(instancePtr, elementSetId, numlay)) return field0Id = ecCreateField(instancePtr) if (.not.(ecSetField1dArray(instancePtr, field0Id, numlay))) return field1Id = ecCreateField(instancePtr) if (.not.(ecSetField1dArray(instancePtr, field1Id, numlay))) return itemId = ecCreateItem(instancePtr) if (.not.ecSetItemRole(instancePtr, itemId, itemType_source)) return if (.not.ecSetItemType(instancePtr, itemId, accessType_fileReader)) return if (.not.ecSetItemQuantity(instancePtr, itemId, quantityId)) return if (.not.ecSetItemElementSet(instancePtr, itemId, elementSetId)) return if (.not.ecSetItemSourceT0Field(instancePtr, itemId, field0Id)) return if (.not.ecSetItemSourceT1Field(instancePtr, itemId, field1Id)) return valueptr => ecSupportFindItem(instancePtr, itemId) ! ===== finish initialization of Fields ===== ! Read the first two records into tEcItem%sourceT0FieldPtr and tEcItem%sourceT1FieldPtr. select case (fileReaderPtr%ofType) case (provFile_t3D) if (.not.ecArcinfoAndCurviReadBlock(fileReaderPtr, fileReaderPtr%fileHandle, valueptr, 0, numlay, 1)) return if (.not.ecArcinfoAndCurviReadBlock(fileReaderPtr, fileReaderPtr%fileHandle, valueptr, 1, numlay, 1)) return case (provFile_bc) if (.not.ecBCReadLine(fileReaderPtr, valueptr%sourceT0FieldPtr%arr1dPtr, valueptr%sourceT0FieldPtr%timesteps)) return if (.not.ecBCReadLine(fileReaderPtr, valueptr%sourceT1FieldPtr%arr1dPtr, valueptr%sourceT1FieldPtr%timesteps)) return case default call setECMessage("ERROR: ec_provider::ecProviderCreatet3DItems: Unknown file type.") return end select ! Add successfully created source Items to the FileReader if (.not.ecAddItemToFileReader(instancePtr, fileReaderPtr%id, valueptr%id)) return success = .true. end function ecProviderCreatet3DItems !> Create subproviders, which create source Items and their contained types. !! meteo1.f90: read1polylin function ecProviderCreatePolyTimItems(instancePtr, fileReaderPtr) result(success) logical :: success !< function status type(tEcInstance), pointer :: instancePtr !< intent(in) type(tEcFileReader), pointer :: fileReaderPtr !< intent(inout) ! real(hp), dimension(:), allocatable :: xs !< x-coordinates of support points real(hp), dimension(:), allocatable :: ys !< y-coordinates of support points real(hp), dimension(:), allocatable :: zs !< z/sigma-coordinates of support points integer, dimension(:), allocatable :: mask !< support point mask array (for polytime ElementSet) integer :: n_points !< number of support points integer :: n_signals !< Number of forcing signals created (at most n_signals==n_points, but warn if n_signals==0) character(len=132) :: rec !< a read line integer :: i, j !< loop counters integer :: istat !< status of read operation integer :: L !< helper index character(len=4) :: tex !< helper string for constructin file names character(len=maxFileNameLen) :: filename !< helper string containing subprovider file name. character(len=maxFileNameLen) :: plipointlbl !< temporary name of current pli-point in bct context character(len=maxFileNameLen) :: polyline_name !< polyline name read from pli-file logical :: exists !< helper boolian, indicating file existence integer :: id !< dummy, catches ids which are not used integer :: k_yyyymmdd !< calculated Gregorian calender date, serving as reference date integer :: quantityId, elementSetId, fieldId, itemId, subconverterId, connectionId, BCBlockID integer :: wind_x, wind_y integer :: magnitude, discharge, waterlevel, slope, crossing, maxLay integer :: latest_uniform_item type(tEcItem), pointer :: itemPT type(tEcItem), pointer :: itemt3D type(tEcItem), pointer :: sourceItem integer, dimension(:), allocatable :: itemIDList integer :: vectormax integer :: lastTimItemID logical :: is_tim, is_cmp, is_tim3d, is_qh type(tEcFileReader), pointer :: fileReaderPtr2 ! ! initialization success = .false. itemPT => null() sourceItem => null() itemt3D => null() maxlay = 0 vectormax = 1 lastTimItemID = -1 ! ! Skip the lead comment lines plus one additional line. do read(fileReaderPtr%fileHandle, '(a)', iostat = istat) rec if (istat /= 0) then call setECMessage("ERROR: ec_provider::ecProviderCreatePolyTimItems: Unexpected end of file.") return end if if (rec(1:1) /= '*') exit end do ! Read the polyline name from the first line of a tekal-block polyline_name = trim(adjustl(rec)) ! Read the number of support points. read(fileReaderPtr%fileHandle, *, iostat = istat) n_points if (istat /= 0) then call setECMessage("ERROR: ec_provider::ecProviderCreatePolyTimItems: Unable to read the number of points.") return end if ! Sanity check if (n_points < 2) then call setECMessage("ERROR: ec_provider::ecProviderCreatePolyTimItems: Less then two support points found.") return end if ! Read the support point coordinate pairs. allocate(xs(n_points)) allocate(ys(n_points)) allocate(mask(n_points)) allocate(itemIDList(n_points)) itemIDList = ec_undef_int mask = 1 do i=1, n_points read(fileReaderPtr%fileHandle, *, iostat = istat) xs(i), ys(i) if (istat /= 0) then call setECMessage("ERROR: ec_provider::ecProviderCreatePolyTimItems: Unable to read a coordinate pair.") return end if enddo ! Construct the poly_tim Item. quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'polytim_item'))) return elementSetId = ecCreateElementSet(instancePtr) if (.not. ecSetElementSetType(instancePtr, elementSetId, elmSetType_cartesian)) return if (.not. ecSetElementSetXArray(instancePtr, elementSetId, xs)) return if (.not. ecSetElementSetYArray(instancePtr, elementSetId, ys)) return if (.not. ecSetElementSetNumberOfCoordinates(instancePtr, elementSetId, n_points)) return fieldId = ecCreateField(instancePtr) itemId = ecCreateItem(instancePtr) if (.not. ecSetItemRole(instancePtr, itemId, itemType_target)) return if (.not. ecSetItemType(instancePtr, itemId, accessType_fileReader)) return if (.not. ecSetItemQuantity(instancePtr, itemId, quantityId)) return if (.not. ecSetItemElementSet(instancePtr, itemId, elementSetId)) return if (.not. ecSetItemTargetField(instancePtr, itemId, fieldId)) return itemPT => ecSupportFindItem(instancePtr, itemId) call jul2ymd(int(fileReaderPtr%tframe%k_refdate + 2400000.5_hp), k_yyyymmdd) ! Init BCBlock for (global) qh-bound is_qh = .false. ! Determine the end of the base of the fileName. L = index(fileReaderPtr%fileName, '.', back = .true.) - 1 ! Create providers at each support point, depending on the availability of specific files. call jul2ymd(int(fileReaderPtr%tframe%k_refdate + 2400000.5_hp), k_yyyymmdd) ! Exceptional case: A single qh-table supplies all support points of the pli-file. filename = fileReaderPtr%fileName(1:L)//'.qh' inquire (file = trim(filename), exist = exists) if (exists) then ! Process a *.qh file. ! Construct a new FileReader id = ecCreateFileReader(instancePtr) if (id == ec_undef_int) then return end if ! Initialize the new FileReader. if (.not. ecProviderInitializeFileReader(instancePtr, id, provFile_qhtable, filename, k_yyyymmdd, fileReaderPtr%tframe%k_timezone)) return ! All fine: is_qh = .true. endif n_signals = 0 ! Record whether at least one child provider is created for this polytim. if (is_qh) then ! Construct a new Converter. subconverterId = ecCreateConverter(instancePtr) ! Determine the source Items. discharge = ecFindItemInFileReader(instancePtr, id, 'discharge') waterlevel = ecFindItemInFileReader(instancePtr, id, 'waterlevel') slope = ecFindItemInFileReader(instancePtr, id, 'slope') crossing = ecFindItemInFileReader(instancePtr, id, 'crossing') if (discharge /= ec_undef_int .and. waterlevel /= ec_undef_int .and. & slope /= ec_undef_int .and. crossing /= ec_undef_int) then !do i=1, n_points ! commented: only one value per polyline ! Initialize the new Converter. if (.not. (ecSetConverterType(instancePtr, subconverterId, convType_qhtable) .and. & ecSetConverterOperand(instancePtr, subconverterId, operand_replace_element) .and. & ecSetConverterInterpolation(instancePtr, subconverterId, interpolate_passthrough) .and. & ecSetConverterElement(instancePtr, subconverterId, 1))) return ! set to 1 from i: only one value per polyline ! Construct a new Connection. connectionId = ecCreateConnection(instancePtr) if (.not. ecSetConnectionConverter(instancePtr, connectionId, subconverterId)) return ! Initialize the new Connection. if (.not. ecAddConnectionSourceItem(instancePtr, connectionId, discharge)) return if (.not. ecAddConnectionSourceItem(instancePtr, connectionId, waterlevel)) return if (.not. ecAddConnectionSourceItem(instancePtr, connectionId, slope)) return if (.not. ecAddConnectionSourceItem(instancePtr, connectionId, crossing)) return if (.not. ecAddConnectionTargetItem(instancePtr, connectionId, itemId)) return if (.not. ecAddItemConnection(instancePtr, itemId, connectionId)) return n_signals = 1 !end do end if else ! .not.is_qh n_signals = 0 do i=1, n_points is_tim = .false. is_cmp = .false. is_tim3d = .false. write(filename,'(a,i4.4,a)') fileReaderPtr%fileName(1:L)//'_',i,'.tim' inquire (file = trim(filename), exist = exists) if (exists) then id = ecCreateFileReader(instancePtr) if (id == ec_undef_int) return if (.not. (ecProviderInitializeFileReader(instancePtr, id, provFile_uniform, filename, k_yyyymmdd, fileReaderPtr%tframe%k_timezone))) return latest_uniform_item = ecFindItemInFileReader(instancePtr, id, 'uniform_item') fileReaderPtr2=>ecSupportFindFileReader(instancePtr, id) is_tim = .true. else write(filename,'(a,i4.4,a)') fileReaderPtr%fileName(1:L)//'_',i,'.cmp' inquire (file = trim(filename), exist = exists) if (exists) then id = ecCreateFileReader(instancePtr) if (id == ec_undef_int) return if (.not. (ecProviderInitializeFileReader(instancePtr, id, provFile_fourier, filename, k_yyyymmdd, fileReaderPtr%tframe%k_timezone))) return fileReaderPtr2=>ecSupportFindFileReader(instancePtr, id) is_cmp = .true. else write(filename,'(a,i4.4,a)') fileReaderPtr%fileName(1:L)//'_',i,'.t3d' inquire (file = trim(filename), exist = exists) if (exists) then id = ecCreateFileReader(instancePtr) if (id == ec_undef_int) return if (.not. (ecProviderInitializeFileReader(instancePtr, id, provFile_t3D, filename, k_yyyymmdd, fileReaderPtr%tframe%k_timezone))) return fileReaderPtr2=>ecSupportFindFileReader(instancePtr, id) is_tim3d = .true. endif ! tim3d-file ? endif ! cmp file ? endif ! tim-file ? if (.not. ecProviderConnectSourceItemsToTargets(instancePtr, is_tim, is_cmp, is_tim3d, id, itemId, i, n_signals, maxlay, itemIDList, lastTimItemID)) then ! ! No sub-FileReader made. mask(i) = 0 endif end do ! loop over support points endif ! switch between qh/cmp/tim if (n_signals <= 0) then call setECMessage("ERROR: ec_provider::ecProviderPolyTimItems: No forcing signals could be attached for polyline file '"//trim(fileReaderPtr%filename)//"'.") return end if if (.not. ecSetElementSetMaskArray(instancePtr, elementSetId, mask)) return ! itemID refers to the source item (providing to the polytim item) for the last support point we came across in the above loop. if (.not. ecProvider3DVectmax(instancePtr, elementSetId, fieldId, mask ,maxlay, n_points, itemIDList, lastTimItemID)) return ! Since the main FileReader's Item is a target, the TimeFrame is not set. ! Add successfully created source Item to the main FileReader if (.not. ecAddItemToFileReader(instancePtr, fileReaderPtr%id, itemPT%id)) return ! ! close pli file close(fileReaderPtr%fileHandle, iostat = istat) success = .true. end function ecProviderCreatePolyTimItems !============================================================================================================== !> Create subproviders, which create source Items and their contained types. !! meteo1.f90: read1polylin function ecProviderCreatePolyTimItemsBC(instancePtr, fileReaderPtr, bctfilename, quantityname) result(success) logical :: success !< function status type(tEcInstance), pointer :: instancePtr !< intent(in) type(tEcFileReader), pointer :: fileReaderPtr !< intent(inout) character(len=*), intent(in) :: bctfilename !< in case of bct-data, we neeed the explicit filename character(len=*), intent(in) :: quantityname !< in case of bct-data, we neeed the explicit quantityname ! real(hp), dimension(:), allocatable :: xs !< x-coordinates of support points real(hp), dimension(:), allocatable :: ys !< y-coordinates of support points real(hp), dimension(:), allocatable :: zs !< z/sigma-coordinates of support points integer, dimension(:), allocatable :: mask !< support point mask array (for polytime ElementSet) integer :: n_points !< number of support points integer :: n_signals !< Number of forcing signals created (at most n_signals==n_points, but warn if n_signals==0) character(len=132) :: rec !< a read line integer :: i, j !< loop counters integer :: istat !< status of read operation integer :: L !< helper index character(len=4) :: tex !< helper string for constructin file names character(len=maxFileNameLen) :: filename !< helper string containing subprovider file name. character(len=maxFileNameLen) :: plipointlbl !< temporary name of current pli-point in bct context character(len=maxFileNameLen) :: polyline_name !< polyline name read from pli-file logical :: exists !< helper boolian, indicating file existence integer :: id !< dummy, catches ids which are not used integer :: k_yyyymmdd !< calculated Gregorian calender date, serving as reference date integer :: quantityId, elementSetId, fieldId, itemId, subconverterId, connectionId, BCBlockID integer :: wind_x, wind_y integer :: magnitude, discharge, waterlevel, slope, crossing, maxLay type(tEcItem), pointer :: itemPT type(tEcItem), pointer :: itemt3D type(tEcItem), pointer :: sourceItem integer, dimension(:), allocatable :: itemIDList integer :: vectormax integer :: lastTimItemID logical :: is_tim, is_cmp, is_tim3d, is_qh type(tEcBCBlock), pointer :: bcBlockPtr type(tEcFileReader), pointer :: fileReaderPtr2 logical :: all_points_are_corr ! ! initialization success = .false. itemPT => null() sourceItem => null() itemt3D => null() maxlay = 0 vectormax = 1 lastTimItemID = -1 ! ! Skip the lead comment lines plus one additional line. do read(fileReaderPtr%fileHandle, '(a)', iostat = istat) rec if (istat /= 0) then call setECMessage("ERROR: ec_provider::ecProviderCreatePolyTimItems: Unexpected end of file.") return end if if (rec(1:1) /= '*') exit end do ! Read the polyline name from the first line of a tekal-block polyline_name = trim(adjustl(rec)) ! Read the number of support points. read(fileReaderPtr%fileHandle, *, iostat = istat) n_points if (istat /= 0) then call setECMessage("ERROR: ec_provider::ecProviderCreatePolyTimItems: Unable to read the number of points.") return end if ! Sanity check if (n_points < 2) then call setECMessage("ERROR: ec_provider::ecProviderCreatePolyTimItems: Less then two support points found.") return end if ! Read the support point coordinate pairs. allocate(xs(n_points)) allocate(ys(n_points)) allocate(mask(n_points)) allocate(itemIDList(n_points)) itemIDList = ec_undef_int mask = 1 do i=1, n_points read(fileReaderPtr%fileHandle, *, iostat = istat) xs(i), ys(i) if (istat /= 0) then call setECMessage("ERROR: ec_provider::ecProviderCreatePolyTimItems: Unable to read a coordinate pair.") return end if enddo call jul2ymd(int(fileReaderPtr%tframe%k_refdate + 2400000.5_hp), k_yyyymmdd) ! Construct the poly_tim Item quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'polytim_item'))) return elementSetId = ecCreateElementSet(instancePtr) if (.not. ecSetElementSetType(instancePtr, elementSetId, elmSetType_cartesian)) return if (.not. ecSetElementSetXArray(instancePtr, elementSetId, xs)) return if (.not. ecSetElementSetYArray(instancePtr, elementSetId, ys)) return if (.not. ecSetElementSetNumberOfCoordinates(instancePtr, elementSetId, n_points)) return fieldId = ecCreateField(instancePtr) itemId = ecCreateItem(instancePtr) if (.not. ecSetItemRole(instancePtr, itemId, itemType_target)) return if (.not. ecSetItemType(instancePtr, itemId, accessType_fileReader)) return if (.not. ecSetItemQuantity(instancePtr, itemId, quantityId)) return if (.not. ecSetItemElementSet(instancePtr, itemId, elementSetId)) return if (.not. ecSetItemTargetField(instancePtr, itemId, fieldId)) return itemPT => ecSupportFindItem(instancePtr, itemId) all_points_are_corr = .true. ! Init BCBlock for (global) qh-bound is_qh = .false. n_signals = 0 ! Record whether at least one child provider is created for this polytim. bcBlockId = ecCreateBCBlock(InstancePtr) bcBlockPtr=>ecSupportFindBCBlock(instancePtr, bcBlockId) write(plipointlbl,'(a)') trim(polyline_name) if (ecProviderInitializeBCBlock(InstancePtr, bcBlockId, k_yyyymmdd, fileReaderPtr%tframe%k_timezone, id, bctfilename, quantityname, plipointlbl, istat)) then is_qh = (bcBlockPtr%func == BC_FUNC_QHTABLE) ! if a polylinename exist as a label without a number ! it might refer to a qh forcing ! SHRL: TODO: check of this is always the case else ! .not.is_qh n_signals = 0 do i=1, n_points is_tim = .false. is_cmp = .false. is_tim3d = .false. ! Process a *.tim file. bcBlockId = ecCreateBCBlock(InstancePtr) bcBlockPtr=>ecSupportFindBCBlock(instancePtr, bcBlockId) ! id van de filereader write(plipointlbl,'(a,i4.4)') trim(polyline_name)//'_', i ! using polyline_name from tekal-block if (.not. ecProviderInitializeBCBlock(InstancePtr, bcBlockId, k_yyyymmdd, fileReaderPtr%tframe%k_timezone, id, bctfilename, quantityname, plipointlbl, istat)) then !call setECMessage("WARNING: ec_provider::ecProviderPolyTimItems: Error initializing EC Block.") mask(i) = 0 mask(i) = 0 cycle endif if (bcBlockPtr%func == BC_FUNC_HARMOCORR .or. bcBlockPtr%func == BC_FUNC_ASTROCORR) then ecAtLeastOnPointIsCorrection = .true. ! n_signals = n_signals + 1 cycle else all_points_are_corr = .false. endif is_tim = (bcBlockPtr%func == BC_FUNC_TSERIES) is_cmp = ((bcBlockPtr%func == BC_FUNC_HARMONIC) .or. (bcBlockPtr%func == BC_FUNC_ASTRO)) is_tim3d = (bcBlockPtr%func == BC_FUNC_TIM3D) if (.not. ecProviderConnectSourceItemsToTargets(instancePtr, is_tim, is_cmp, is_tim3d, id, itemId, i, n_signals, maxlay, itemIDList, lastTimItemID)) then ! ! No sub-FileReader made. mask(i) = 0 endif end do ! loop over support points endif ! switch between qh/cmp/tim if (ecAtLeastOnPointIsCorrection) then if (all_points_are_corr) then success = .true. else ! TODO: error: bc file currently should only contain corrections endif return endif if (n_signals <= 0) then call setECMessage("ERROR: ec_provider::ecProviderPolyTimItemsBC: No forcing signals could be attached for polyline file '"//trim(fileReaderPtr%filename)//"'.") return end if if (is_qh) then ! Construct a new Converter. subconverterId = ecCreateConverter(instancePtr) ! Determine the source Items. discharge = ecFindItemInFileReader(instancePtr, id, 'discharge') waterlevel = ecFindItemInFileReader(instancePtr, id, 'waterlevel') slope = ecFindItemInFileReader(instancePtr, id, 'slope') crossing = ecFindItemInFileReader(instancePtr, id, 'crossing') if (discharge /= ec_undef_int .and. waterlevel /= ec_undef_int .and. & slope /= ec_undef_int .and. crossing /= ec_undef_int) then !do i=1, n_points ! commented: only one value per polyline ! Initialize the new Converter. if (.not. (ecSetConverterType(instancePtr, subconverterId, convType_qhtable) .and. & ecSetConverterOperand(instancePtr, subconverterId, operand_replace_element) .and. & ecSetConverterInterpolation(instancePtr, subconverterId, interpolate_passthrough) .and. & ecSetConverterElement(instancePtr, subconverterId, 1))) return ! set to 1 from i: only one value per polyline ! Construct a new Connection. connectionId = ecCreateConnection(instancePtr) if (.not. ecSetConnectionConverter(instancePtr, connectionId, subconverterId)) return ! Initialize the new Connection. if (.not. ecAddConnectionSourceItem(instancePtr, connectionId, discharge)) return if (.not. ecAddConnectionSourceItem(instancePtr, connectionId, waterlevel)) return if (.not. ecAddConnectionSourceItem(instancePtr, connectionId, slope)) return if (.not. ecAddConnectionSourceItem(instancePtr, connectionId, crossing)) return if (.not. ecAddConnectionTargetItem(instancePtr, connectionId, itemId)) return if (.not. ecAddItemConnection(instancePtr, itemId, connectionId)) return n_signals = 1 !end do end if endif if (.not. ecSetElementSetMaskArray(instancePtr, elementSetId, mask)) return ! itemID refers to the source item (providing to the polytim item) for the last support point we came across in the above loop. if (.not. ecProvider3DVectmax(instancePtr, elementSetId, fieldId, mask ,maxlay, n_points, itemIDList, lastTimItemID)) return ! Since the main FileReader's Item is a target, the TimeFrame is not set. ! Add successfully created source Item to the main FileReader if (.not. ecAddItemToFileReader(instancePtr, fileReaderPtr%id, itemPT%id)) return ! ! close pli file close(fileReaderPtr%fileHandle, iostat = istat) success = .true. end function ecProviderCreatePolyTimItemsBC !============================================================================================================== function ecProvider3DVectmax(instancePtr, elementSetId, fieldId, mask ,maxlay, n_points, itemIDList, sourceItemID) result (success) implicit none logical :: success type(tEcInstance), pointer :: instancePtr !< intent(in) integer, intent(in) :: elementSetId integer, intent(in) :: fieldId integer, intent(in) :: mask(:) integer, intent(in) :: maxlay integer, intent(in) :: n_points integer, dimension(:), allocatable :: itemIDList integer, intent(in) :: sourceItemID real(hp), dimension(:), allocatable :: zs !< z/sigma-coordinates of support points integer :: i, magnitude, vectormax type(tEcItem), pointer :: itemt3D, sourceItem success = ecSetElementSetMaskArray(instancePtr, elementSetId, mask) !! Dit wordt een aparte routine: additonal 3D stuff if ( maxlay > 0 ) then ! t3D type allocate(zs(maxlay*n_points)) zs = ec_undef_hp do i=1, n_points magnitude = itemIDList(i) if ( magnitude /= ec_undef_int ) then Itemt3D => ecSupportFindItem(instancePtr, magnitude) zs((i-1)*maxlay+1:(i-1)*maxlay+itemt3D%elementSetPtr%nCoordinates) = itemt3D%elementSetPtr%x end if end do if (success) success = ecSetElementSetZArray(instancePtr, elementSetId, zs) if (.not. ecSetField1dArray(instancePtr, fieldId, n_points*maxlay)) then success = .false. end if else ! Determine vectormax. sourceItem => ecSupportFindItem(instancePtr, sourceItemID) if ( associated(sourceItem) ) then if (sourceItem%elementSetPtr%ofType == elmSetType_scalar) then vectormax = sourceItem%elementSetPtr%nCoordinates end if end if ! Allocate the target Field. if (.not. ecSetField1dArray(instancePtr, fieldId, n_points*vectormax)) then success = .false. end if end if end function ecProvider3DVectmax !============================================================================================================== function ecProviderConnectSourceItemsToTargets(instancePtr, is_tim, is_cmp, is_tim3d, fileReaderId, targetItemId, targetIndex, n_signals, maxlay, itemIDList, lastTimItemID) result(itemFound) logical :: itemFound type(tEcInstance), pointer :: instancePtr !< intent(in) integer :: fileReaderId ! file reader id integer :: targetItemId ! target item id integer :: targetIndex ! index in target item's values integer :: n_signals, maxlay ! INOUT integer, dimension(:) :: itemIDList ! INOUT integer, intent(inout) :: lastTimItemID logical :: is_tim, is_cmp, is_tim3d integer :: subconverterId, magnitude, j, connectionId, nr_fourier_items, anItemId type(tEcItem), pointer :: itemt3D itemFound = .false. if (is_tim) then ! Construct a new Converter. subconverterId = ecCreateConverter(instancePtr) ! Determine the source Items. magnitude = ecFindItemInFileReader(instancePtr, fileReaderId, 'uniform_item') lastTimItemID = magnitude if (magnitude /= ec_undef_int) then ! Initialize the new Converter. if (.not. (ecSetConverterType(instancePtr, subconverterId, convType_uniform) .and. & ecSetConverterOperand(instancePtr, subconverterId, operand_replace_element) .and. & ecSetConverterInterpolation(instancePtr, subconverterId, interpolate_timespace) .and. & ecSetConverterElement(instancePtr, subconverterId, targetIndex))) return ! Construct a new Connection. connectionId = ecCreateConnection(instancePtr) if (.not. ecSetConnectionConverter(instancePtr, connectionId, subconverterId)) return ! Initialize the new Connection. if (.not. ecAddConnectionSourceItem(instancePtr, connectionId, magnitude)) return if (.not. ecAddConnectionTargetItem(instancePtr, connectionId, targetItemId)) return if (.not. ecAddItemConnection(instancePtr, targetItemId, connectionId)) return n_signals = n_signals + 1 else ! error handling end if itemFound = .true. end if if (is_cmp) then ! Construct a new Converter. subconverterId = ecCreateConverter(instancePtr) ! Determine the source Items. nr_fourier_items = ecGetFileReaderNumberOfItems(instancePtr, fileReaderId) ! Initialize the new Converter. if (.not. (ecSetConverterType(instancePtr, subconverterId, convType_fourier) .and. & ecSetConverterOperand(instancePtr, subconverterId, operand_replace_element) .and. & ecSetConverterInterpolation(instancePtr, subconverterId, interpolate_passthrough) .and. & ecSetConverterElement(instancePtr, subconverterId, targetIndex))) return ! Construct a new Connection. connectionId = ecCreateConnection(instancePtr) if (.not. ecSetConnectionConverter(instancePtr, connectionId, subconverterId)) return ! Initialize the new Connection. do j=1, nr_fourier_items anItemId = ecFindXthItemInFileReader(instancePtr, fileReaderId, j) if (.not. ecAddConnectionSourceItem(instancePtr, connectionId, anItemId)) return end do if (.not. ecAddConnectionTargetItem(instancePtr, connectionId, targetItemId)) return if (.not. ecAddItemConnection(instancePtr, targetItemId, connectionId)) return n_signals = n_signals + 1 itemFound = .true. end if if (is_tim3d) then ! Construct a new Converter. subconverterId = ecCreateConverter(instancePtr) ! Determine the source Items. magnitude = ecFindItemInFileReader(instancePtr, fileReaderId, 'quant') ! update maximum number of layers Itemt3D => ecSupportFindItem(instancePtr, magnitude) maxlay = max(maxlay,Itemt3D%elementSetPtr%nCoordinates) ItemIDList(targetIndex) = magnitude ! Initialize the new Converter. if (.not. (ecSetConverterType(instancePtr, subconverterId, convType_uniform) .and. & ecSetConverterOperand(instancePtr, subconverterId, operand_replace_element) .and. & ecSetConverterInterpolation(instancePtr, subconverterId, interpolate_time) .and. & ecSetConverterElement(instancePtr, subconverterId, targetIndex))) return ! Construct a new Connection. connectionId = ecCreateConnection(instancePtr) if (.not. ecSetConnectionConverter(instancePtr, connectionId, subconverterId)) return ! Initialize the new Connection. if (.not. ecAddConnectionSourceItem(instancePtr, connectionId, magnitude)) return if (.not. ecAddConnectionTargetItem(instancePtr, connectionId, targetItemId)) return if (.not. ecAddItemConnection(instancePtr, targetItemId, connectionId)) return n_signals = n_signals + 1 itemFound = .true. endif end function ecProviderConnectSourceItemsToTargets ! ======================================================================= !> Create source Items and their contained types, based on a spiderweb file's header. !! meteo1.f90: reaspwheader function ecProviderCreateSpiderwebItems(instancePtr, fileReaderPtr) result(success) logical :: success !< function status type(tEcInstance), pointer :: instancePtr !< intent(in) type(tEcFileReader), pointer :: fileReaderPtr !< intent(inout) ! integer :: quantityId !< helper variable integer :: elementSetId !< helper variable integer :: field0Id !< helper variable integer :: field1Id !< helper variable integer :: itemId !< helper variable integer :: n_cols !< helper variable integer :: n_rows !< helper variable real(hp) :: missingValue !< helper variable real(hp) :: radius !< helper variable character(len=maxNameLen) :: radius_unit !< helper variable type(tEcItem), pointer :: item1 !< Item containing quantity1 type(tEcItem), pointer :: item2 !< Item containing quantity2 type(tEcItem), pointer :: item3 !< Item containing quantity3 character(len=maxNameLen) :: rec !< helper variable ! success = .false. item1 => null() item2 => null() item3 => null() ! rec = ecSpiderwebAndCurviFindInFile(fileReaderPtr%fileHandle, 'n_cols') read(rec, *) n_cols n_cols = n_cols + 1 ! an additional column for 360 == 0 degrees rec = ecSpiderwebAndCurviFindInFile(fileReaderPtr%fileHandle, 'n_rows') read(rec, *) n_rows n_rows = n_rows + 1 ! an additional row for the eye rec = ecSpiderwebAndCurviFindInFile(fileReaderPtr%fileHandle, 'NODATA_value') read(rec, *) missingValue rec = ecSpiderwebAndCurviFindInFile(fileReaderPtr%fileHandle, 'spw_radius') read(rec, *) radius rec = ecSpiderwebAndCurviFindInFile(fileReaderPtr%fileHandle, 'spw_rad_unit') read(rec, *) radius_unit ! ! One common ElementSet. elementSetId = ecCreateElementSet(instancePtr) if (.not. (ecSetElementSetType(instancePtr, elementSetId, elmSetType_spw) .and. & ecSetElementSetRadius(instancePtr, elementSetId, radius, radius_unit) .and. & ecSetElementSetRowsCols(instancePtr, elementSetId, n_rows, n_cols))) then success = .false. end if ! ! ===== quantity1: wind_speed ===== quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'windspeed') .and. & ecSetQuantityUnits(instancePtr, quantityId, trim(ecSpiderwebAndCurviFindInFile(fileReaderPtr%fileHandle, 'unit1'))))) then success = .false. end if field0Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field0Id, n_cols*n_rows) .and. & ecSetFieldMissingValue(instancePtr, field0Id, missingValue))) then success = .false. end if field1Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field1Id, n_cols*n_rows) .and. & ecSetFieldMissingValue(instancePtr, field1Id, missingValue))) then success = .false. end if itemId = ecCreateItem(instancePtr) if (.not. ( ecSetItemRole(instancePtr, itemId, itemType_source) .and. & ecSetItemType(instancePtr, itemId, accessType_fileReader) .and. & ecSetItemQuantity(instancePtr, itemId, quantityId) .and. & ecSetItemElementSet(instancePtr, itemId, elementSetId) .and. & ecSetItemSourceT0Field(instancePtr, itemId, field0Id) .and. & ecSetItemSourceT1Field(instancePtr, itemId, field1Id))) then success = .false. else item1 => ecSupportFindItem(instancePtr, itemId) end if ! ===== quantity2: wind_from_direction ===== quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'winddirection') .and. & ecSetQuantityUnits(instancePtr, quantityId, trim(ecSpiderwebAndCurviFindInFile(fileReaderPtr%fileHandle, 'unit2'))))) then success = .false. end if field0Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field0Id, n_cols*n_rows) .and. & ecSetFieldMissingValue(instancePtr, field0Id, missingValue))) then success = .false. end if field1Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field1Id, n_cols*n_rows) .and. & ecSetFieldMissingValue(instancePtr, field1Id, missingValue))) then success = .false. end if itemId = ecCreateItem(instancePtr) if (.not. ( ecSetItemRole(instancePtr, itemId, itemType_source) .and. & ecSetItemType(instancePtr, itemId, accessType_fileReader) .and. & ecSetItemQuantity(instancePtr, itemId, quantityId) .and. & ecSetItemElementSet(instancePtr, itemId, elementSetId) .and. & ecSetItemSourceT0Field(instancePtr, itemId, field0Id) .and. & ecSetItemSourceT1Field(instancePtr, itemId, field1Id))) then success = .false. else item2 => ecSupportFindItem(instancePtr, itemId) end if ! ===== quantity3: p_drop ===== quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'p_drop') .and. & ecSetQuantityUnits(instancePtr, quantityId, trim(ecSpiderwebAndCurviFindInFile(fileReaderPtr%fileHandle, 'unit3'))))) then success = .false. end if field0Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field0Id, n_cols*n_rows) .and. & ecSetFieldMissingValue(instancePtr, field0Id, missingValue))) then success = .false. end if field1Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field1Id, n_cols*n_rows) .and. & ecSetFieldMissingValue(instancePtr, field1Id, missingValue))) then success = .false. end if itemId = ecCreateItem(instancePtr) if (.not. ( ecSetItemRole(instancePtr, itemId, itemType_source) .and. & ecSetItemType(instancePtr, itemId, accessType_fileReader) .and. & ecSetItemQuantity(instancePtr, itemId, quantityId) .and. & ecSetItemElementSet(instancePtr, itemId, elementSetId) .and. & ecSetItemSourceT0Field(instancePtr, itemId, field0Id) .and. & ecSetItemSourceT1Field(instancePtr, itemId, field1Id))) then success = .false. else item3 => ecSupportFindItem(instancePtr, itemId) end if ! ===== finish initialization of Fields ===== ! Read the first two records into tEcItem%sourceT0FieldPtr and tEcItem%sourceT1FieldPtr. rewind(unit=fileReaderPtr%fileHandle) success = ecSpiderAndCurviAndArcinfoReadToBody(fileReaderPtr%fileHandle) if (success) then success = ecSpiderwebReadBlock(fileReaderPtr, item1, item2, item3, 0, n_cols, n_rows) end if if (success) then success = ecSpiderwebReadBlock(fileReaderPtr, item1, item2, item3, 1, n_cols, n_rows) end if ! Add successfully created source Items to the FileReader if (success) success = ecAddItemToFileReader(instancePtr, fileReaderPtr%id, item1%id) if (success) success = ecAddItemToFileReader(instancePtr, fileReaderPtr%id, item2%id) if (success) success = ecAddItemToFileReader(instancePtr, fileReaderPtr%id, item3%id) end function ecProviderCreateSpiderwebItems ! ======================================================================= !> Create source Items and their contained types, based on a NetCDF file's header. function ecProviderCreateNetcdfItems(instancePtr, fileReaderPtr, quantityName) result(success) logical :: success !< function status type(tEcInstance), pointer :: instancePtr !< intent(in) type(tEcFileReader), pointer :: fileReaderPtr !< intent(inout) character(maxNameLen), intent(in) :: quantityName !< name of quantity to read ! integer :: ierror !< return value of NetCDF function calls integer :: idvar_rainfall !< id as obtained from NetCDF integer :: idvar_coord !< id as obtained from NetCDF integer :: idvar_time !< id as obtained from NetCDF integer :: ndims !< helper variable, containing the number of dimensions integer, dimension(:), allocatable :: dimids !< ids of a variable's dimensions integer, dimension(:), allocatable :: coord_ids !< helper variable character(NF90_MAX_NAME), dimension(:), allocatable :: names !< helper variable, containing dimension names integer :: i !< loop counter integer :: fgd_id !< var_id for elementset X or latitude integer :: sgd_id !< var_id for elementset Y or longitude integer :: fgd_grid_type !< helper variable for consistency check on grid_type integer :: sgd_grid_type !< helper variable for consistency check on grid_type integer :: grid_type !< elmSetType enum integer :: fgd_size !< number of grid points in first dimension integer :: sgd_size !< number of grid points in second dimension real(hp), dimension(:,:), allocatable :: fgd_data !< coordinate data along first dimension's axis real(hp), dimension(:), allocatable :: fgd_data_1d !< coordinate data along first dimension's axis real(hp), dimension(:,:), allocatable :: sgd_data !< coordinate data along second dimension's axis real(hp), dimension(:), allocatable :: sgd_data_1d !< coordinate data along second dimension's axis real(hp) :: fdg_miss !< missing data value in first dimension real(hp) :: sdg_miss !< missing data value in second dimension character(NF90_MAX_NAME) :: units !< helper variable for variable's units character(NF90_MAX_NAME) :: coord_name !< helper variable character(NF90_MAX_NAME), dimension(:), allocatable :: coord_names !< helper variable character(NF90_MAX_NAME) :: name !< helper variable integer :: quantityId !< helper variable integer :: elementSetId !< helper variable integer :: field0Id !< helper variable integer :: field1Id !< helper variable integer :: itemId !< helper variable integer :: istat !< helper variable type(tEcItem), pointer :: itemPtr !< Item containing quantity integer :: t0t1 !< t0 / t1 switch logical :: dummy !< temp integer, dimension(:), allocatable :: dim_sizes !< helper variable ! success = .false. itemPtr => null() fgd_grid_type = ec_undef_int sgd_grid_type = ec_undef_int grid_type = ec_undef_int units = '' coord_name = '' name = '' ndims = 0 ! ============================================================================= ! Find the Quantity corresponding to quantityName. (configurable in the future) ! ============================================================================= ! TODO: Check on standard variable name. ! TODO: Check on variable's long_name. ! Check for variable named "rainfall". FEWS calls this "Rainfall" if (.not. ecSupportNetcdfCheckError(nf90_inq_varid(fileReaderPtr%fileHandle, 'Rainfall', idvar_rainfall), "inq_varid Rainfall", fileReaderPtr%filename)) return ! ============================================= ! Inspect for the presence of grid information. ! ============================================= ! Find the string of whitespace spearated data variables in the "coordinates" attribute of the variable "Rainfall". ! i.e. "lat lon time" ierror = nf90_get_att(fileReaderPtr%fileHandle, idvar_rainfall, "coordinates", coord_name) ! Find sizes of the dimensions of the data variables in the "coordinates" attribute of the variable "Rainfall". ierror = nf90_inquire_variable(fileReaderPtr%fileHandle, idvar_rainfall, ndims=ndims) allocate(dimids(ndims)) ierror = nf90_inquire_variable(fileReaderPtr%fileHandle, idvar_rainfall, dimids=dimids) allocate(dim_sizes(ndims)) if (ndims > 1) then do i=1, ndims ierror = nf90_inquire_dimension(fileReaderPtr%fileHandle, dimids(i), len=dim_sizes(i)) end do end if ! Identify the grid type and couple the data variables to their sizes. allocate(coord_names(ndims)) coord_names = '' read(coord_name, *,iostat=istat) ( coord_names(i), i=1,ndims ) ! As temporary leniency, we will tolerate refrainment of time coordinate specification. if (istat .ne. 0) then coord_names(ndims) = 'time' endif allocate(coord_ids(ndims)) if (ndims > 1) then do i=1, ndims ierror = nf90_inq_varid(fileReaderPtr%fileHandle, coord_names(i), coord_ids(i)) name = '' ! NetCDF fails to overwrite the entire string, so re-initialize each iteration. ierror = nf90_get_att(fileReaderPtr%fileHandle, coord_ids(i), "standard_name", name) if (strcmpi(name, 'projection_x_coordinate')) then fgd_id = coord_ids(i) ierror = nf90_inquire_dimension(fileReaderPtr%fileHandle, dimids(i), len=fgd_size) fgd_grid_type = elmSetType_cartesian end if if (strcmpi(name, 'projection_y_coordinate')) then sgd_id = coord_ids(i) ierror = nf90_inquire_dimension(fileReaderPtr%fileHandle, dimids(i), len=sgd_size) sgd_grid_type = elmSetType_cartesian end if if (strcmpi(name, 'longitude')) then fgd_id = coord_ids(i) ierror = nf90_inquire_dimension(fileReaderPtr%fileHandle, dimids(i), len=sgd_size) fgd_grid_type = elmSetType_spheric end if if (strcmpi(name, 'latitude')) then sgd_id = coord_ids(i) ierror = nf90_inquire_dimension(fileReaderPtr%fileHandle, dimids(i), len=fgd_size) sgd_grid_type = elmSetType_spheric end if end do if (fgd_grid_type /= ec_undef_int .and. sgd_grid_type /= ec_undef_int .and. fgd_grid_type == sgd_grid_type) then grid_type = fgd_grid_type else call setECMessage("ERROR: ec_provider::ecProviderCreateNetcdfItems: Inconsistent grid type between dimensions.") end if end if ! ===================== ! Create the ElementSet ! ===================== elementSetId = ecCreateElementSet(instancePtr) if (grid_type == ec_undef_int) then dummy = ecSetElementSetNumberOfCoordinates(instancePtr, elementSetId, 0) else allocate(fgd_data(fgd_size,sgd_size), sgd_data(fgd_size,sgd_size)) ierror = nf90_get_var(fileReaderPtr%fileHandle, fgd_id, fgd_data, start=(/1,1/), count=(/fgd_size,sgd_size/)) ierror = nf90_get_var(fileReaderPtr%fileHandle, sgd_id, sgd_data, start=(/1,1/), count=(/fgd_size,sgd_size/)) allocate(fgd_data_1d(fgd_size*sgd_size), sgd_data_1d(fgd_size*sgd_size)) fgd_data_1d = reshape(fgd_data, (/fgd_size*sgd_size/)) sgd_data_1d = reshape(sgd_data, (/fgd_size*sgd_size/)) if (grid_type == elmSetType_cartesian) then if (.not. (ecSetElementSetType(instancePtr, elementSetId, grid_type) .and. & ecSetElementSetXArray(instancePtr, elementSetId, fgd_data_1d) .and. & ecSetElementSetYArray(instancePtr, elementSetId, sgd_data_1d) .and. & ecSetElementSetRowsCols(instancePtr, elementSetId, sgd_size, fgd_size) .and. & ecSetElementSetNumberOfCoordinates(instancePtr, elementSetId, fgd_size*sgd_size))) then return end if else if (grid_type == elmSetType_spheric) then if (.not. (ecSetElementSetType(instancePtr, elementSetId, grid_type) .and. & ecSetElementSetLatitudeArray(instancePtr, elementSetId, fgd_data_1d) .and. & ecSetElementSetLongitudeArray(instancePtr, elementSetId, sgd_data_1d) .and. & ecSetElementSetRowsCols(instancePtr, elementSetId, sgd_size, fgd_size) .and. & ecSetElementSetNumberOfCoordinates(instancePtr, elementSetId, fgd_size*sgd_size))) then return end if end if end if ! =================== ! Create the Quantity ! =================== quantityId = ecCreateQuantity(instancePtr) ierror = nf90_get_att(fileReaderPtr%fileHandle, idvar_rainfall, 'units', units) if (.not. (ecSetQuantityName(instancePtr, quantityId, 'Rainfall') .and. & ! ReadBlock uses this name! ecSetQuantityUnits(instancePtr, quantityId, units))) then return end if ! ======================== ! Create the source Fields ! ======================== ! --- Determine missingDataValue --- if (.not. ecSupportNetcdfCheckError(nf90_get_att(fileReaderPtr%fileHandle, fgd_id, "_FillValue", fdg_miss), "reading _FillValue", fileReaderPtr%fileName)) then if (.not. ecSupportNetcdfCheckError(nf90_get_att(fileReaderPtr%fileHandle, NF90_GLOBAL, "NF90_FILL_DOUBLE", fdg_miss), "reading _FillValue", fileReaderPtr%fileName)) then fdg_miss = ec_undef_hp end if end if if (.not. ecSupportNetcdfCheckError(nf90_get_att(fileReaderPtr%fileHandle, fgd_id, "_FillValue", sdg_miss), "reading _FillValue", fileReaderPtr%fileName)) then if (.not. ecSupportNetcdfCheckError(nf90_get_att(fileReaderPtr%fileHandle, NF90_GLOBAL, "NF90_FILL_DOUBLE", sdg_miss), "reading _FillValue", fileReaderPtr%fileName)) then sdg_miss = ec_undef_hp end if end if ! field0Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field0Id, fgd_size*sgd_size) .and. & ecSetFieldMissingValue(instancePtr, field0Id, fdg_miss))) then return end if field1Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field1Id, fgd_size*sgd_size) .and. & ecSetFieldMissingValue(instancePtr, field1Id, sdg_miss))) then return end if ! ================== ! Create source Item ! ================== itemId = ecCreateItem(instancePtr) if (.not. ( ecSetItemRole(instancePtr, itemId, itemType_source) .and. & ecSetItemType(instancePtr, itemId, accessType_fileReader) .and. & ecSetItemQuantity(instancePtr, itemId, quantityId) .and. & ecSetItemElementSet(instancePtr, itemId, elementSetId) .and. & ecSetItemSourceT0Field(instancePtr, itemId, field0Id) .and. & ecSetItemSourceT1Field(instancePtr, itemId, field1Id))) then return else itemPtr => ecSupportFindItem(instancePtr, itemId) end if ! ===== finish initialization of Fields ===== ! Read the first two records into tEcItem%sourceT0FieldPtr and tEcItem%sourceT1FieldPtr. ! Add successfully created source Item to the FileReader success = ecNetcdfReadNextBlock(fileReaderPtr, itemPtr, 0) if (success) then ! trick: set fieldT1%timesteps to fieldT0%timesteps, so that T1 will be filled with the block after T0's block. itemPtr%sourceT1FieldPtr%timesteps = itemPtr%sourceT0FieldPtr%timesteps success = ecNetcdfReadNextBlock(fileReaderPtr, itemPtr, 1) end if ! Add successfully created source Items to the FileReader if (success) success = ecAddItemToFileReader(instancePtr, fileReaderPtr%id, itemPtr%id) end function ecProviderCreateNetcdfItems ! ======================================================================= !> Create source Items and their contained types, based on NetCDF file header. function ecProviderCreateWaveNetcdfItems(instancePtr, fileReaderPtr, quantityName) result(success) logical :: success !< function status type(tEcInstance), pointer :: instancePtr !< intent(in) type(tEcFileReader), pointer :: fileReaderPtr !< intent(inout) character(maxNameLen), intent(in) :: quantityName !< name of quantity to read ! integer :: i !< loop counter integer :: ierror !< return value of function calls integer :: iddim_netelem !< id as obtained from NetCDF integer :: idvar_x !< id as obtained from NetCDF integer :: idvar_q !< id as obtained from NetCDF integer :: n !< number of values integer :: quantityId !< helper variable integer :: elementSetId !< helper variable integer :: field0Id !< helper variable integer :: field1Id !< helper variable integer :: itemId !< helper variable integer :: t0t1 !< indicates whether the 0 or the 1 field is read. -1: choose yourself logical :: local_success !< when the return flag should not be influenced real(hp) :: dmiss !< missing data value type(tEcItem), pointer :: item character(20) :: name !< character(NF90_MAX_NAME) :: string !< read from NetCDF file character(300) :: message ! ! body success = .false. item => null() dmiss = -999.0_hp elementSetId = ecCreateElementSet(instancePtr) ! ! Cartesian or spheric ! WARNING: elementSetType must be set before elementSetNumberOfCoordinates (why???) ! string = ' ' ierror = nf90_get_att(fileReaderPtr%fileHandle, nf90_global, 'grid_mapping', string); success = ecSupportNetcdfCheckError(ierror, "get_att global grid_mapping", fileReaderPtr%filename) if (string == 'wgs84') then success = ecSetElementSetType(instancePtr, elementSetId, elmSetType_spheric) else success = ecSetElementSetType(instancePtr, elementSetId, elmSetType_cartesian) endif ! ! Read the number of values on the NetCDF file; this should exactly match the number of cells in dflowfm ! ierror = nf90_inq_dimid(fileReaderPtr%fileHandle, 'nFlowElemWithBnd', iddim_netelem); success = ecSupportNetcdfCheckError(ierror, "inq_dimid netelem", fileReaderPtr%fileName) ierror = nf90_inquire_dimension(fileReaderPtr%fileHandle, iddim_netelem, string, n); success = ecSupportNetcdfCheckError(ierror, "inq_dim netelem", fileReaderPtr%fileName) success = ecElementSetSetNumberOfCoordinates(instancePtr, elementSetId, n) ! ! Use quantityName ! string = ' ' ierror = nf90_inq_varid(fileReaderPtr%fileHandle, quantityName, idvar_q); success = ecSupportNetcdfCheckError(ierror, "inq_varid " // quantityName, fileReaderPtr%filename) ierror = nf90_get_att(fileReaderPtr%fileHandle, idvar_q, 'units', string); success = ecSupportNetcdfCheckError(ierror, "inq_att " // quantityName, fileReaderPtr%filename) if (.not.success) return ! quantityId = ecCreateQuantity(instancePtr) if (.not. (ecSetQuantityName(instancePtr, quantityId, quantityName) .and. & ecSetQuantityUnits(instancePtr, quantityId, trim(string)))) then success = .false. endif field0Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field0Id, n) .and. & ecSetFieldMissingValue(instancePtr, field0Id, dmiss))) then success = .false. endif field1Id = ecCreateField(instancePtr) if (.not. (ecSetField1dArray(instancePtr, field1Id, n) .and. & ecSetFieldMissingValue(instancePtr, field1Id, dmiss))) then success = .false. endif itemId = ecCreateItem(instancePtr) if (.not. ( ecSetItemRole(instancePtr, itemId, itemType_source) .and. & ecSetItemType(instancePtr, itemId, accessType_fileReader) .and. & ecSetItemQuantity(instancePtr, itemId, quantityId) .and. & ecSetItemElementSet(instancePtr, itemId, elementSetId) .and. & ecSetItemSourceT0Field(instancePtr, itemId, field0Id) .and. & ecSetItemSourceT1Field(instancePtr, itemId, field1Id))) then success = .false. else item => ecSupportFindItem(instancePtr, itemId) endif ! ! finish initialization of Fields ! Read the first two records into tEcItem%sourceT0FieldPtr and tEcItem%sourceT1FieldPtr ! if (success) then t0t1 = 0 success = ecNetcdfReadBlock(fileReaderPtr, item, t0t1, n) endif if (success) then t0t1 = 1 local_success = ecNetcdfReadBlock(fileReaderPtr, item, t0t1, n) endif ! Add successfully created source Item to the FileReader if (success) then success = ecAddItemToFileReader(instancePtr, fileReaderPtr%id, item%id) endif if (.not. success) then write(message,'(2a)') "ERROR: ec_provider::ecProviderCreateWaveNetcdfItems: unable to create sourceItem for ", quantityName call setECMessage(trim(message)) endif end function ecProviderCreateWaveNetcdfItems ! ======================================================================= !> Set the TimeFrame's properties, based on a file header when present, or to refdat and hardcoded values otherwise. function ecProviderInitializeTimeFrame(fileReaderPtr, refdat, tzone) result(success) logical :: success !< function status type(tEcFileReader), pointer :: fileReaderPtr !< intent(inout) integer, intent(in) :: refdat !< Kernel's reference date, format: Gregorian yyyymmdd real(hp), intent(in) :: tzone !< Kernel's timezone. ! success = .false. ! fileReaderPtr%tframe%k_refdate = dble(ymd2jul(refdat) - 2400000.5_hp) fileReaderPtr%tframe%k_timezone = tzone fileReaderPtr%tframe%ec_refdate = fileReaderPtr%tframe%k_refdate ! select case(fileReaderPtr%ofType) case (provFile_undefined) call setECMessage("ERROR: ec_provider::ecProviderInitializeTimeFrame: No file type specified.") case (provFile_uniform, provFile_unimagdir) success = ecUniInitializeTimeFrame(fileReaderPtr) case (provFile_svwp) call setECMessage("ERROR: ec_provider::ecProviderInitializeTimeFrame: Unsupported file type.") case (provFile_svwp_weight) call setECMessage("ERROR: ec_provider::ecProviderInitializeTimeFrame: Unsupported file type.") case (provFile_arcinfo) success = ecDefaultInitializeTimeFrame(fileReaderPtr) case (provFile_spiderweb) success = ecDefaultInitializeTimeFrame(fileReaderPtr) case (provFile_curvi) success = ecDefaultInitializeTimeFrame(fileReaderPtr) case (provFile_curvi_weight) call setECMessage("ERROR: ec_provider::ecProviderInitializeTimeFrame: Unsupported file type.") case (provFile_triangulation) call setECMessage("ERROR: ec_provider::ecProviderInitializeTimeFrame: Unsupported file type.") case (provFile_triangulationmagdir) call setECMessage("ERROR: ec_provider::ecProviderInitializeTimeFrame: Unsupported file type.") case (provFile_qhtable) fileReaderPtr%tframe%ec_refdate = fileReaderPtr%tframe%k_refdate ! The nr_timesteps is infinity. success = .true. case (provFile_poly_tim) ! TODO : For now the top-level fileReader (*.pli) needs no TimeFrame. fileReaderPtr%tframe%ec_refdate = fileReaderPtr%tframe%k_refdate ! TODO : The total number of available timesteps is not yet known here. success = .true. case (provFile_fourier) ! Boundary condition components file containing: either period[minute] or component name, amplitude[m], phase[degree]. ! Time_unit explicit in the Converter, as values at a moment in time are calculated from the seeds in the file, not read from file as records. fileReaderPtr%tframe%ec_refdate = fileReaderPtr%tframe%k_refdate ! The nr_timesteps is infinity, as values at a moment in time are calculated from the seeds in the file, not read from file as records. success = .true. case (provFile_grib) call setECMessage("ERROR: ec_provider::ecProviderInitializeTimeFrame: Unsupported file type.") case (provFile_netcdf) success = ecNetcdfInitializeTimeFrame(fileReaderPtr) case (provFile_t3D) success = ecDefaultInitializeTimeFrame(fileReaderPtr) case (provFile_bc) ! Filereader was created thru a BC-instance. This instance has a timeunit (netcdf-style) property if (fileReaderPtr%bc%func == BC_FUNC_TSERIES) then success = ecSupportTimestringToUnitAndRefdate(fileReaderPtr%bc%timeunit, & fileReaderPtr%tframe%ec_timestep_unit, fileReaderPtr%tframe%ec_refdate) else success = .true. endif case default call setECMessage("ERROR: ec_provider::ecProviderInitializeTimeFrame: Unknown file type.") end select end function ecProviderInitializeTimeFrame ! ======================================================================= !> Set the TimeFrame's properties, based on a default time string format. !! meteo1.f90: reaspwtim function ecDefaultInitializeTimeFrame(fileReaderPtr) result(success) logical :: success !< function status type(tEcFileReader), pointer :: fileReaderPtr !< intent(in) ! character(len=maxNameLen) :: keyword !< helper variable character(len=maxNameLen) :: rec !< helper variable character(len=maxNameLen) :: prev_rec !< helper variable integer :: indx !< helper variable integer :: year !< helper variable integer :: month !< helper variable integer :: day !< helper variable real(hp) :: time_steps !< helper variable ! success = .false. keyword = 'TIME' prev_rec = ' ' ! rewind(unit=fileReaderPtr%fileHandle) ! rec = ecFindInFile(fileReaderPtr%fileHandle, keyword) if (len_trim(rec) > 0) then ! Determine the timestep unit and reference date for the time data in the file. if (.not. ecSupportTimestringToUnitAndRefdate(rec, fileReaderPtr%tframe%ec_timestep_unit, fileReaderPtr%tframe%ec_refdate)) return else call setECMessage("ERROR: ec_provider::ecDefaultInitializeTimeFrame: Unable to identify the first data block.") return end if ! ========================================= ! Determine the total number of time steps. ! ========================================= !!! This read action is too slow on multi-gigabyte files. Commented for now. !!!do !!! prev_rec = rec !!! rec = ecFindInFile(fileReaderPtr%fileHandle, keyword) !!! if (len_trim(rec) == 0) exit !!!end do !!!if (len_trim(prev_rec) /= 0) then !!! ! Read and convert the timesteps to seconds. !!! if (.not. ecGetTimesteps(prev_rec, time_steps, .false.)) return !!! ! Store the total number of available timesteps in the ecFileReader's ecTimeFrame. !!! fileReaderPtr%tframe%nr_timesteps = time_steps !!!else !!! call setECMessage("ERROR: ec_provider::ecDefaultInitializeTimeFrame: Unable to read the stop time.") !!!end if success = .true. end function ecDefaultInitializeTimeFrame ! ======================================================================= !> Set the TimeFrame's properties, based on a uni* file's header. !! Headerless column-based ASCII file. !! uniform: time[minute], wind_x[m s-1](, wind_y[m s-1]) !! unimagdir: time[minute], windspeed[m s-1], winddirection[degree] function ecUniInitializeTimeFrame(fileReaderPtr) result(success) logical :: success !< function status type(tEcFileReader), pointer :: fileReaderPtr !< FileReader to initialize ! real(hp) :: time_steps !< number of timesteps real(hp) :: prev_time_steps !< helper variable for determining tstop ! success = .false. fileReaderPtr%tframe%ec_refdate = fileReaderPtr%tframe%k_refdate ! Obtain the total number of timesteps in this file. rewind(unit=fileReaderPtr%fileHandle) !!! This read action is too slow on multi-gigabyte files. Commented for now. !!!time_steps = 0.0_hp !!!do !!! prev_time_steps = time_steps !!! if (.not. ecUniReadTimeSteps(fileReaderPtr, time_steps)) exit ! TODO: EB: do we really need to count all time steps beforehand? !!!end do !!!! Convert timesteps from minutes to seconds. !!!fileReaderPtr%tframe%nr_timesteps = prev_time_steps * 60.0_hp success = .true. ! ecSupportTimestringToUnitAndRefdate(units, fileReaderPtr%tframe%ec_timestep_unit, fileReaderPtr%tframe%ec_refdate) end function ecUniInitializeTimeFrame ! ======================================================================= !> Set the TimeFrame's properties, based on a NetCDF file. function ecNetcdfInitializeTimeFrame(fileReaderPtr) result(success) use netcdf ! logical :: success !< function status type(tEcFileReader), pointer :: fileReaderPtr !< intent(in) ! integer :: nVariables !< number of variables in NetCDF file integer :: time_id !< integer id of variable with standard_name "time" integer :: i !< loop counter character(len=NF90_MAX_NAME) :: name !< name of a variable character(len=NF90_MAX_NAME) :: units !< units attribute of a variable integer, dimension(1) :: dimid !< integer id of time variable's dimension variable integer :: length !< number of time steps integer :: istat !< status of allocation operation ! success = .false. nVariables = 0 time_id = ec_undef_int ! ! Determine the total number of variables inside the NetCDF file. if (.not. ecSupportNetcdfCheckError(nf90_inquire(fileReaderPtr%fileHandle, nVariables=nVariables), "obtain nVariables", fileReaderPtr%fileName)) return ! ! Inspect the standard_name attribute of all variables to find "time" and store that variable's id. do i=1, nVariables name = '' ! NetCDF does not completely overwrite a string, so re-initialize. if (.not. ecSupportNetcdfCheckError(nf90_get_att(fileReaderPtr%fileHandle, i, "standard_name", name), "obtain standard_name", fileReaderPtr%fileName)) then cycle endif if (strcmpi(name, 'time')) then time_id = i exit end if end do if (time_id == ec_undef_int) then call setECMessage("ERROR: ec_provider::ecNetcdfInitializeTimeFrame: Unable to find variable with standard_name: time.") return end if ! ! Determine the timestep unit and reference date for the time data in the NetCDF file. ! Surprisingly, the reference date is part of the "units" attribute. units = '' ! NetCDF does not completely overwrite a string, so re-initialize. if (.not. ecSupportNetcdfCheckError(nf90_get_att(fileReaderPtr%fileHandle, time_id, "units", units), "obtain units", fileReaderPtr%fileName)) return if (.not. ecSupportTimestringToUnitAndRefdate(units, fileReaderPtr%tframe%ec_timestep_unit, fileReaderPtr%tframe%ec_refdate)) return ! ! Determine the total number of timesteps. if (.not. ecSupportNetcdfCheckError(nf90_inquire_variable(fileReaderPtr%fileHandle, time_id, dimids=dimid), "obtain time dimension ids", fileReaderPtr%fileName)) return if (.not. ecSupportNetcdfCheckError(nf90_inquire_dimension(fileReaderPtr%fileHandle, dimid(1), len=length), "obtain time dimension length", fileReaderPtr%fileName)) return fileReaderPtr%tframe%nr_timesteps = 0.0_hp + length allocate(fileReaderPtr%tframe%times(length), stat = istat) ! Store the times at which data is available. if (.not. ecSupportNetcdfCheckError(nf90_get_var(fileReaderPtr%fileHandle, time_id, fileReaderPtr%tframe%times, start=(/1/), count=(/length/)), "obtain time data", fileReaderPtr%fileName)) return ! success = .true. end function ecNetcdfInitializeTimeFrame ! ======================================================================= !> Read the header of a Arcinfo file. !! Does not check whether all expected information is present in the header. !! meteo1:: readarcinfoheader function readarcinfoheader(minp, mmax, nmax, x0, y0, dxa, dya, dmiss) result(success) logical :: success !< function status integer :: minp !< integer, intent(out) :: mmax !< integer, intent(out) :: nmax !< real(hp), intent(out) :: x0 !< real(hp), intent(out) :: y0 !< real(hp), intent(out) :: dxa !< real(hp), intent(out) :: dya !< real(hp), intent(out) :: dmiss !< ! integer :: jacornerx integer :: jacornery integer :: k integer :: l integer :: equal_sign_index character(len=maxNameLen) :: rec ! success = .false. jacornerx = 0 jacornery = 0 ! rewind(unit = minp) 10 continue read (minp, '(A)', end = 100) rec if (index(rec, '### START OF HEADER') > 0) then ! new d3dflow header 20 continue read (minp, '(A)', end = 100) rec equal_sign_index = index(rec, '=') L = equal_sign_index + 1 ! if (index(rec, 'NODATA_value') > 0 .and. index(rec, 'NODATA_value') < equal_sign_index) then read (rec(L:), *, err = 106) dmiss endif ! if (index(rec, 'n_cols') > 0 .and. index(rec, 'n_cols') < equal_sign_index) then read (rec(L:), *, err = 101) mmax endif ! if (index(rec, 'n_rows') > 0 .and. index(rec, 'n_rows') < equal_sign_index) then read (rec(L:), *, err = 102) nmax endif ! if (index(rec, 'x_llcenter') > 0 .and. index(rec, 'x_llcenter') < equal_sign_index) then read (rec(L:), *, err = 103) x0 endif ! if (index(rec, 'x_llcorner') > 0 .and. index(rec, 'x_llcorner') < equal_sign_index) then read (rec(L:), *, err = 103) x0 jacornerx = 1 endif ! if (index(rec, 'y_llcenter') > 0 .and. index(rec, 'y_llcenter') < equal_sign_index) then read (rec(L:), *, err = 104) y0 endif ! if (index(rec, 'y_llcorner') > 0 .and. index(rec, 'y_llcorner') < equal_sign_index) then read (rec(L:), *, err = 104) y0 jacornery = 1 endif ! if (index(rec, 'dx') > 0 .and. index(rec, 'dx') < equal_sign_index) then read (rec(L:), *, err = 105) dxa endif ! if (index(rec, 'dy') > 0 .and. index(rec, 'dy') < equal_sign_index) then read (rec(L:), *, err = 105) dya endif ! if (index(rec, '### END OF HEADER') == 0) then ! new d3dflow header goto 20 endif else ! if (rec(1:1)=='*' .or. rec(2:2)=='*') goto 10 read (rec(13:), *, err = 101) mmax read (minp, '(A)', end = 100) rec read (rec(13:), *, err = 102) nmax ! read (minp, '(A)', end = 100) rec read (rec(13:), *, err = 103) x0 jacornerx = 0 if (index(rec, 'corner')/=0) jacornerx = 1 ! read (minp, '(A)', end = 100) rec read (rec(13:), *, err = 104) y0 jacornery = 0 if (index(rec, 'corner')/=0) jacornery = 1 ! read (minp, '(A)', end = 100) rec l = index(rec, 'cellsize') + 8 k = numbersonline(rec(l:)) if (k==1) then read (rec(13:), *, err = 105) dxa dya = dxa jacornery = jacornerx else read (rec(13:), *, err = 105) dxa, dya endif read (minp, '(A)', end = 100) rec read (rec(13:), *, err = 106) dmiss ! ! Data in an arcinfo grid file is always defined in the cell centres. ! Data is assumed to be given at the points x0+i*dx,y0+j*dy ! If the x0/y0 line contains the word corner, the corner coordinates ! have been specified in the file. Therefore, shift x0 and y0 by half ! a grid cell to the cell centres. ! endif ! if (jacornerx == 1) x0 = x0 + dxa/2 if (jacornery == 1) y0 = y0 + dya/2 ! success = .true. return ! ! error handling ! 100 continue call setECMessage('ERROR: Unexpected end of file while reading header of arcinfo file') goto 999 101 continue call setECMessage('ERROR: Looking for ncols (arc-info), but getting: ',trim(rec)) goto 999 102 continue call setECMessage('ERROR: Looking for nrows (arc-info), but getting: ',trim(rec)) goto 999 103 continue call setECMessage('ERROR: Looking for xll (arc-info), but getting: ',trim(rec)) goto 999 104 continue call setECMessage('ERROR: Looking for yll (arc-info), but getting: ',trim(rec)) goto 999 105 continue call setECMessage('ERROR: Looking for cellsize (dx, dy) (arc-info), but getting: ',trim(rec)) goto 999 106 continue call setECMessage('ERROR: Looking for missing value (arc-info), but getting: ',trim(rec)) goto 999 999 continue success = .false. end function readarcinfoheader ! ======================================================================= function numbersonline(rec) integer :: numbersonline character(*) :: rec ! integer :: i integer :: istarti integer :: leeg integer :: lend ! numbersonline = 0 leeg = 1 lend = len_trim(rec) do i = 1, lend if (index(rec(i:i), ' ')==0) then ! there is nothing here if (leeg==1) then leeg = 0 istarti = i numbersonline = numbersonline + 1 endif else leeg = 1 endif enddo end function numbersonline end module m_ec_provider