module m_ec_bcreader use precision use m_ec_parameters use m_ec_stringbuffer use m_ec_support use m_ec_message use m_ec_typedefs use m_alloc use multi_file_io use string_module implicit none private public :: ecInitBC public :: ecBCReadLine public :: ecBCBlockCreate public :: jakeyvalue contains !subroutine uppr(r) !implicit none !character(len=*), intent(inout) :: r !integer :: i, ch !do i=1,len_trim(r) ! ch = ichar(r(i:i)) ! if(ch>=97) then ! if(ch<=122) then ! r(i:i)=achar(ch-97+65) ! endif ! endif !enddo !return !end subroutine uppr logical function ecInitBC (fname, fhandle, quantityName, plilabel, bc, iostat) result (success) implicit none ! Open a .bc file with external forcings ! Scrolls the headers until it finds quantity, location and .... are matched ! Advances the filepointer upto the first line of data, setting all necesary fields of bc-object ! Leaves the file open and return the filehandle fhandle (in a bc-object) to caller character(len=*), intent(in) :: fname integer(kind=8), intent(in) :: fhandle character(len=*), intent(in) :: quantityName character(len=*), intent(in) :: plilabel type (tEcBCBlock), intent(inout) :: bc integer, optional, intent(out) :: iostat type (tEcFileReader), pointer :: fileReaderPtr success = .false. bc%qname = quantityName bc%bcname = plilabel call str_upper(bc%qname,len(trim(bc%qname))) call str_upper(bc%bcname,len(trim(bc%bcname))) ! if (scanbc(bc, fhandle, iostat)) then ! parsing the open bc-file bc%fname = fname ! store original filename for later referece (?) bc%fhandle = fhandle allocate(bc%columns(bc%numcols)) else call mf_close(fhandle) return endif success=.true. end function ecInitBC logical function scanbc(bc, fhandle, iostat) result (found) implicit none type (tEcBCBlock), intent(inout) :: bc integer (kind=8), intent(in) :: fhandle integer, intent(out) :: iostat character*(255) :: rec integer :: reclen integer :: commentpos character*(1000) :: keyvaluestr ! all key-value pairs in one header integer :: posfs integer :: nfld integer :: nq logical :: jablock, jaheader integer :: reallocstat integer :: lineno integer (kind=8) :: savepos integer :: iostatloc found = .false. iostat = EC_UNKNOWN_ERROR lineno = 0 savepos = 0 do if (mf_eof(fhandle)) then ! print *,'End of BC FILE' ! reached the end of the bc-file, but combination was not found iostat = EC_EOF return endif ! read(fhandle,'(a200)', iostat=iostat) rec call mf_read(fhandle,rec,savepos) iostatloc = 0 ! mf_read always ok? if (iostatloc /= 0) then iostat = iostatloc return ! beter break? endif lineno = lineno + 1 reclen = len_trim(rec) ! deal with various comment delimiters commentpos = index(rec,'//') if (commentpos>0) reclen = min(reclen,commentpos-1) commentpos = index(rec,'%') if (commentpos>0) reclen = min(reclen,commentpos-1) commentpos = index(rec,'#') if (commentpos>0) reclen = min(reclen,commentpos-1) commentpos = index(rec,'*') if (commentpos>0) reclen = min(reclen,commentpos-1) commentpos = index(rec,'!') if (commentpos>0) reclen = min(reclen,commentpos-1) if (len_trim(rec(1:reclen))>0) then ! skip empty lines if (index(rec,'[forcing]')>0) then ! new boundary chapter jaheader = .true. ! switching to header mode keyvaluestr = ',' jablock=.false. nfld = 0 ! count the number of fields in this header block nq = 0 ! count the (maximum) number of quantities in this block else if (jaheader) then posfs = index(rec(1:reclen),'=') ! key value pair ? if (posfs>0) then call replace_char(rec,9,32) ! replace tabs by spaces, header key-value pairs only nfld = nfld + 1 ! count the number of lines in the header file ! Create a lengthy string of ',key,value,key,value.....' call str_upper(rec(1:posfs-1)) ! all keywords uppercase , not case sensitive if (index(rec(1:posfs-1),'QUANTITY')>0) then nq = nq + 1 endif ! keyvaluestr = trim(keyvaluestr)//(trim(adjustl(rec(1:posfs-1))))//',' & ! //(trim(adjustl(rec(posfs+1:reclen))))//',' keyvaluestr = trim(keyvaluestr)//''''// (trim(adjustl(rec(1:posfs-1))))//''',''' & //(trim(adjustl(rec(posfs+1:reclen))))//''',' else ! switch to datamode call str_upper(keyvaluestr,len(trim(keyvaluestr))) ! case insensitive format if (jakeyvalue(keyvaluestr,'NAME',trim(bc%bcname))) then if (jakeyvalue(keyvaluestr,'FUNCTION','ASTRONOMIC').or.jakeyvalue(keyvaluestr,'FUNCTION','ASTRONOMIC-CORRECTION')) then ! Check for harmonic or astronomic components jablock = .true. jablock = jablock .and. jakeyvalue(keyvaluestr,'QUANTITY','ASTRONOMIC COMPONENT') jablock = jablock .and. jakeyvalue(keyvaluestr,'QUANTITY',trim(bc%qname)//' '//'AMPLITUDE') jablock = jablock .and. jakeyvalue(keyvaluestr,'QUANTITY',trim(bc%qname)//' '//'PHASE') elseif (jakeyvalue(keyvaluestr,'FUNCTION','HARMONIC').or.jakeyvalue(keyvaluestr,'FUNCTION','HARMONIC-CORRECTION')) then ! Check for harmonic or harmonic components jablock = .true. jablock = jablock .and. jakeyvalue(keyvaluestr,'QUANTITY','HARMONIC COMPONENT') jablock = jablock .and. jakeyvalue(keyvaluestr,'QUANTITY',trim(bc%qname)//' '//'AMPLITUDE') jablock = jablock .and. jakeyvalue(keyvaluestr,'QUANTITY',trim(bc%qname)//' '//'PHASE') elseif (jakeyvalue(keyvaluestr,'FUNCTION','QH')) then ! Check for qh jablock = .true. jablock = jablock .and. jakeyvalue(keyvaluestr,'QUANTITY',trim(bc%qname)//' '//'WATERLEVEL') jablock = jablock .and. jakeyvalue(keyvaluestr,'QUANTITY',trim(bc%qname)//' '//'DISCHARGE') elseif (jakeyvalue(keyvaluestr,'FUNCTION','T3D')) then ! Check for timeseries on sigma- or z-levels jablock = jakeyvalue(keyvaluestr,'QUANTITY',bc%qname) elseif (jakeyvalue(keyvaluestr,'FUNCTION','TIMESERIES')) then ! Check for timeseries jablock = jakeyvalue(keyvaluestr,'QUANTITY',bc%qname) endif if (jablock) then ! block confirmed call processhdr(bc,nfld,nq,keyvaluestr) ! dumb translation of bc-object metadata call checkhdr(bc) ! check on the contents of the bc-object call mf_backspace(fhandle, savepos) ! Rewind the first line with data found = .true. iostat = EC_NOERR return else ! location was found, but not all required meta data was present iostat = EC_METADATA_INVALID endif ! Right quantity else ! location was found, but data was missing/invalid iostat = EC_DATA_NOTFOUND endif ! Right label jaheader = .false. ! No, we are NOT reading a header if (jablock) then ! Yes, we are in the bc-block of interest endif ! found matching block endif ! switch to datamode endif ! in header mode (data lines are ignored) endif ! not a new '[forcing]' item endif ! non-empty string enddo ! read/scan loop end function scanbc function jakeyvalue(keyvaluestr,key,value) result (jafound) implicit none logical :: jafound character(len=*), intent(in) :: keyvaluestr character(len=*), intent(in) :: key character(len=*), intent(in) :: value jafound = (index(keyvaluestr,','''//trim(key)//''','''//trim(value)//''',')>0) end function jakeyvalue !> Given a character string of key-value pairs gathered from a header block, !> extract all relevant fields (and find the block of requested quantity) !> and fill a bc-object subroutine processhdr(bc,nfld,nq,keyvaluestr) implicit none type (tEcBCBlock), intent(inout) :: bc !< boundary condition instance integer, intent(in) :: nfld !< number of fields in the header integer, intent(in) :: nq !< number of quantities (quantity block) character(len=*), intent(in) :: keyvaluestr !< string containing header fields as key-value pairs integer :: ifld integer :: i integer :: iostat character(len=60), allocatable :: hdrkeys(:) !< All keys from header character(len=60), allocatable :: hdrvals(:) !< All values from header character(len=60) :: dumstr integer :: ipos, npos integer :: iq, iq_sel if (allocated(hdrkeys)) then deallocate(hdrkeys) endif if (allocated(hdrvals)) then deallocate(hdrvals) endif if (allocated(bc%quantity%jacolumn)) then deallocate(bc%quantity%jacolumn) endif if (allocated(bc%quantity%vect2col)) then deallocate(bc%quantity%vect2col) endif allocate(hdrkeys(nfld),hdrvals(nfld)) allocate(bc%quantity%jacolumn(nq)) allocate(bc%quantity%vect2col(nq)) bc%quantity%vect2col = -1 hdrvals='' hdrkeys='' read(keyvaluestr,*,iostat=iostat) dumstr,(hdrkeys(ifld),hdrvals(ifld),ifld=1,nfld) iq = 0 iq_sel = 0 do ifld=1,nfld select case (trim(adjustl(hdrkeys(ifld)))) case ('QUANTITY') iq = iq + 1 ! count quantities, corresponds with column numbers [iq] if (trim(hdrvals(ifld))==bc%qname) then ! detected quantity of interest bc%quantity%jacolumn(iq)=.true. iq_sel = iq_sel + 1 bc%quantity%vect2col(iq_sel) = iq else bc%quantity%jacolumn(iq)=.false. select case (bc%func) case (BC_FUNC_TSERIES, BC_FUNC_TIM3D) if (trim(hdrvals(ifld))=='TIME') then ! special check on the time field bc%timecolumn = iq endif case (BC_FUNC_HARMONIC, BC_FUNC_ASTRO, BC_FUNC_HARMOCORR, BC_FUNC_ASTROCORR, BC_FUNC_CMP3D) if (trim(hdrvals(ifld))=='HARMONIC COMPONENT') then ! harmonic component bc%quantity%astro_component_column = iq endif if (trim(hdrvals(ifld))=='ASTRONOMIC COMPONENT') then ! astronomic component label bc%quantity%astro_component_column = iq endif if (trim(hdrvals(ifld))==trim(bc%qname)//' AMPLITUDE') then ! amplitude field for astronomic/harmonic components bc%quantity%astro_amplitude_column = iq endif if (trim(hdrvals(ifld))==trim(bc%qname)//' PHASE') then ! phase field for astronomic/harmonic components bc%quantity%astro_phase_column = iq endif case (BC_FUNC_QHTABLE) if (trim(hdrvals(ifld))==trim(bc%qname)//' WATERLEVEL') then ! waterlevel field for qh-boundary bc%quantity%qh_waterlevel_column = iq endif if (trim(hdrvals(ifld))==trim(bc%qname)//' DISCHARGE') then ! discharge field for qh-boundary bc%quantity%qh_discharge_column = iq endif end select ! functions endif case ('UNIT') if (bc%quantity%jacolumn(iq)) then bc%quantity%unit = trim(hdrvals(ifld)) endif if (iq==bc%timecolumn) then ! Is this the unit of time ? bc%timeunit = trim(hdrvals(ifld)) ! store timeunit string in this bc instance endif case ('FUNCTION') select case (trim(adjustl(hdrvals(ifld)))) case ('TIMESERIES') bc%func = BC_FUNC_TSERIES case ('HARMONIC') bc%func = BC_FUNC_HARMONIC case ('ASTRONOMIC') bc%func = BC_FUNC_ASTRO case ('HARMONIC-CORRECTION') bc%func = BC_FUNC_HARMOCORR case ('ASTRONOMIC-CORRECTION') bc%func = BC_FUNC_ASTROCORR case ('T3D') bc%func = BC_FUNC_TIM3D case ('C3D') bc%func = BC_FUNC_CMP3D case ('QHTABLE') bc%func = BC_FUNC_QHTABLE end select case ('OFFSET') ! QUESTION: Offset and Factor are Unit-related right? if (bc%quantity%jacolumn(iq)) then read(hdrvals(ifld),*) bc%quantity%offset endif case ('FACTOR') if (bc%quantity%jacolumn(iq)) then read(hdrvals(ifld),*) bc%quantity%factor endif case ('VERTICAL POSITION') ! if (bc%quantity%jacolumn(iq)) then ! if multiple vertical levels are permitted ! npos=0 ! if (len_trim(hdrvals(ifld))>0) then ! npos = count([(verify(hdrvals(ifld)(i:i),', ')>0 & ! .and.verify(hdrvals(ifld)(i-1:i-1),', ')==0, i=2,len_trim(hdrvals(ifld)))]) + 1 ! endif ! allocate(bc%quantity%vertndx(npos)) ! read(hdrvals(ifld),*) (bc%quantity%vertndx(ipos),ipos=1,npos) ! endif read(hdrvals(ifld),*) bc%quantity%vertndx case ('VERTICAL POSITION SPECIFICATION') npos=0 if (len_trim(hdrvals(ifld))>0) then npos = count([(verify(hdrvals(ifld)(i:i),', ')>0 & .and.verify(hdrvals(ifld)(i-1:i-1),', ')==0, i=2,len_trim(hdrvals(ifld)))]) + 1 endif allocate(bc%vp(npos)) read(hdrvals(ifld),*) (bc%vp(ipos),ipos=1,npos) ! globally store ALL vertical positions case ('MISSING VALUE DEFINITION') read(hdrvals(ifld),*) bc%missing case ('TIME-INTERPOLATION') select case (trim(adjustl(hdrvals(ifld)))) case ('LINEAR') bc%vertint = BC_TIMEINT_LIN case ('BLOCK-TO') bc%vertint = BC_TIMEINT_BTO case ('BLOCK-FROM') bc%vertint = BC_TIMEINT_BFROM end select case ('VERTICAL INTERPOLATION') select case (trim(adjustl(hdrvals(ifld)))) case ('LINEAR') bc%vertint = BC_VERTINT_LIN case ('LOG') bc%vertint = BC_VERTINT_LOG case ('BLOCK') bc%vertint = BC_VERTINT_BLOCK end select case ('VERTICAL POSITION TYPE') IF (index(hdrvals(ifld),'PERCEN')+index(hdrvals(ifld),'BED')>0) then hdrvals(ifld)='PERCBED' endif select case (trim(adjustl(hdrvals(ifld)))) case ('PERCBED') bc%vptyp = BC_VPTYP_PERCBED case ('ZDATUM') bc%vptyp = BC_VPTYP_ZDATUM case ('BEDSURF') bc%vptyp = BC_VPTYP_BEDSURF case ('PERCSURF') bc%vptyp = BC_VPTYP_PERCSURF case ('ZBED') bc%vptyp = BC_VPTYP_ZBED case ('ZSURF') bc%vptyp = BC_VPTYP_ZSURF end select end select enddo bc%numcols = iq deallocate(hdrkeys) deallocate(hdrvals) end subroutine processhdr !> Given a filled bc-object, scrutinize its content for completeness, validity and consistency subroutine checkhdr(bc) implicit none type (tEcBCBlock), intent(in) :: bc ! ! end subroutine checkhdr !> Read the next record from a *.bc file. !> Requests a line from the EcBC object's stringbuffer block, advancing its pointer in the block function ecBCReadLine(fileReaderPtr, values, time_steps, recout, eof) result(success) implicit none logical :: success !< function status type(tEcFileReader), pointer :: fileReaderPtr real(hp), optional, intent(inout) :: time_steps !< number of time steps of duration: seconds real(hp), dimension(:), optional, intent(inout) :: values !< vector of values for a BC_FUNC_TSERIES character(len=*), optional, intent(out) :: recout !< line prepared from input for caller logical, optional, intent(out) :: eof !< reading failed, but only because eof ! type(tEcBCBlock), pointer :: bcPtr integer :: n_col !< number of columns in the file, inferred from the number of quantity blocks in the header integer :: n_col_time !< position of the dedicated time-column in the line, deduced from the header character(256) :: rec !< content of a line integer :: istat !< status of read operation integer :: i, j !< loop counters integer :: reclen integer :: commentpos integer(kind=8):: savepos !< saved position in file, for mf_read to enabled rewinding real(kind=hp) :: ec_timesteps ! to read in source time from file block bcPtr => fileReaderPtr%bc success = .false. n_col = bcPtr%numcols ! Number of quantities = columns specified in the header n_col_time = bcPtr%timecolumn ! Rank number of the colunm presumably holding the time (not necessarily the first) if(present(recout)) then recout = '' ! initialize return string to empty endif rec = '' reclen = len_trim(rec) if (present(eof)) then eof = .false. endif do while(reclen==0) if (mf_eof(bcPtr%fhandle)) then call setECMessage("INFO: ec_bcreader::ecBCReadLine: Datablock end (eof) has been reached of: "//trim(bcPtr%fname)) if (present(eof)) then eof = .true. endif return endif call mf_read(bcPtr%fhandle,rec,savepos) reclen = len_trim(rec) ! deal with various comment delimiters, RL: skip it if performance is affected ! commentpos = index(rec,'//') ! Scan for various types of comments if (commentpos>0) reclen = min(reclen,commentpos-1) commentpos = index(rec,'%') if (commentpos>0) reclen = min(reclen,commentpos-1) commentpos = index(rec,'#') if (commentpos>0) reclen = min(reclen,commentpos-1) commentpos = index(rec,'*') if (commentpos>0) reclen = min(reclen,commentpos-1) commentpos = index(rec,'!') if (commentpos>0) reclen = min(reclen,commentpos-1) reclen = len_trim(rec(1:reclen)) ! Finally remove trailing spaces if (index(rec,'[forcing]')>0) then call setECMessage("INFO: ec_bcreader::ecBCReadBlock: Datablock end (new [forcing] block) has been reached of: "//trim(bcPtr%fname)) eof = .true. return endif enddo ! read(rec(1:reclen), *, IOSTAT = istat) (values(i), i=1,n_col_time-1), time_steps, (values(i), i=n_col_time, n_col-1) read(rec(1:reclen), *, IOSTAT = istat) BCPtr%columns(1:n_col) if (istat /= 0) then ! error handling, report column number i, field content columns(i) and record rec .... ! TODO: hookup MessageHandlign and print rec and column stats here directly call setECMessage("ERROR: ec_bcreader::ecBCReadBlock: Read failure in file: "//trim(bcPtr%fname)) success = .false. return endif select case (BCPtr%func) case (BC_FUNC_TSERIES, BC_FUNC_TIM3D) read (BCPtr%columns(n_col_time), *) ec_timesteps ! Convert source time to kernel time: time_steps = ecSupportThisTimeToTimesteps(fileReaderPtr%tframe, ec_timesteps) j=0 if (count(BCPtr%quantity%vect2col>0)==0) then ! vect2col is not used to map the read colums, added to the vector in reading order do i=1,n_col if (BCPtr%quantity%jacolumn(i)) then j=j+1 read(BCPtr%columns(i),*,iostat=istat) values(j) if (istat/=0) then ! error handling, report column number i, field content columns(i) and record rec .... call setECMessage("ERROR: ec_bcreader::ecBCReadBlock: Read failure in file: "//trim(bcPtr%fname)) return endif endif enddo else ! vect2col is used to map vector elements to columns (useful if a permutation of columns is desired) do i=1,n_col j=j+1 if (BCPtr%quantity%vect2col(i)>0) then read(BCPtr%columns(BCPtr%quantity%vect2col(i)),*,iostat=istat) values(j) endif enddo endif case (BC_FUNC_ASTRO,BC_FUNC_HARMONIC,BC_FUNC_ASTROCORR,BC_FUNC_HARMOCORR) ! Produce a record of component, amplitude, phase extracted from rec recout = trim(BCPtr%columns(BCPtr%quantity%astro_component_column)) recout = trim(recout)//' '//trim(BCPtr%columns(BCPtr%quantity%astro_amplitude_column)) recout = trim(recout)//' '//trim(BCPtr%columns(BCPtr%quantity%astro_phase_column)) case (BC_FUNC_QHTABLE) ! Produce a record of waterlevel, discharge extracted from rec recout = trim(BCPtr%columns(BCPtr%quantity%qh_discharge_column)) recout = trim(recout)//' '//trim(BCPtr%columns(BCPtr%quantity%qh_waterlevel_column)) end select success = .true. end function ecBCReadLine ! ======================================================================= ! Realloc sub and oldfil function to be removed lateron, ! only included for standalone-testing ! realloc can be found in m_alloc (deltares common) ! oldfil to be replaced by ! ecSupportOpenExistingFile(fhandle, fname) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - integer function oldfil(fname, maxunit, iostat) result(unit) implicit none character(len=*), intent(in) :: fname integer, intent(in) :: maxunit integer, intent(out) :: iostat logical :: opened iostat = 0 ! ran out of free file handles do unit=1, maxunit inquire(unit,opened=opened) if (.not.opened) then exit endif enddo if (unit<=maxunit) then open(unit, file=trim(fname), status='old', iostat=iostat) else iostat = -66 ! ran out of free file handles endif ! ..... return end function oldfil ! ======================================================================= !> Construct a new FileReader with the specified id. !! Failure is indicated by returning a null pointer. function ecBCBlockCreate(bcBlockId) result(bcBlockPtr) type(tEcBCBlock), pointer :: bcBlockPtr !< the new BCBlock, intent(out) integer, intent(in) :: bcBlockId !< unique BCBlock id ! integer :: istat !< allocate() status integer :: i !< loop counter logical :: success !< helper variable ! success = .false. ! ! allocation allocate(bcBlockPtr, stat = istat) if (istat == 0) then ! see what's allocatable in this object and add it here ! ...... end if if (istat /= 0) then call setECMessage("ERROR: ec_bcreader::ecBCBlockCreate: Unable to allocate additional memory.") bcBlockPtr => null() return end if ! initialization bcBlockPtr%id = bcBlockId end function ecBCBlockCreate end module m_ec_bcreader