module input use utilities, only : & open_file, & log_debug, & log_info, & log_warn, & log_error, & dp, & check_range, & check_int_range, & path_max_length implicit none ! Path to config file set via commandline argument character(len=path_max_length) :: config_file_path = "" integer :: & ! Incremental line reading size in characters, must fit max. column name length read_line_buffer_size = 128 integer, parameter :: SEEK_SET = 0, SEEK_CUR = 1, SEEK_END = 2 ! Input prefix character(len=path_max_length), save :: input_dir = "input/" ! Counter for offset in discharge.csv column, as leap year logic etc is ! in outer loops integer :: discharge_offset = 1 namelist / INPUT_PARAMETERS / & input_dir #ifdef with_netcdf ! INPUT FILES AND VARIABLE NAMES ! variables read from ncinfo file character(len=path_max_length) :: NC_GRID = "subbasin_climate_grid.csv" character(len=path_max_length) :: NC_FNAMES(6) ! netCDF file names for 6 variables character(len=path_max_length) :: NC_VNAMES(6) ! netCDF variable names for 6 variables character(len = 100) :: NC_LAT_VNAME = "lat" character(len = 100) :: NC_LON_VNAME = "lon" character(len = 100) :: NC_TIME_VNAME = "time" ! reference year of integer time variable (uninterpreted if < 0) integer :: NC_REF_YEAR = - 1 ! number of days to offset irregardless of NC_REF_YEAR integer :: NC_OFFSET_DAYS = 0 logical :: NC_DEBUG = .true. namelist / NC_CLIMATE_PARAMETERS / & NC_FNAMES, & NC_VNAMES, & NC_LAT_VNAME, & NC_LON_VNAME, & NC_TIME_VNAME, & NC_REF_YEAR, & NC_OFFSET_DAYS, & NC_GRID, & NC_DEBUG ! Attributes character(len=14) :: nc_attr_names(4) = (/'scale_factor ', 'add_offset ', & '_FillValue ', 'missing_value'/) ! all need to be the same length real(dp) :: nc_attr_vals(6, 4) = 0 ! IDs and counters integer :: nc_ids(6), nc_var_ids(6) integer :: nc_nlons, nc_nlats, nc_nrecs integer :: nc_xmin = 10000000 integer :: nc_xmax = 0 integer :: nc_ymin = 10000000 integer :: nc_ymax = 0 integer :: nc_nx, nc_ny ! current record number, ie. day counter from begining of file (reset in input_nc_find_time) integer :: nc_nday = 1 ! number of rows in ncinfo file, set later integer :: nc_nrows = 0 ! max cells per subbasin integer :: nc_mxc = 0 ! allocated with nc_mxc in readinfo ! weights read from input file real(dp), dimension(:, :), allocatable :: nc_weight_sb ! lon/lat read from input file real(dp), dimension(:, :), allocatable :: nc_lon_sb, nc_lat_sb ! indecies of lon / lat in nc files integer, dimension(:, :), allocatable :: nc_x_sb, nc_y_sb ! allocatable data arrays ! sum of weights per subbasin real(dp), dimension(:), allocatable :: nc_wsum_sb ! number of cells per subbasin integer, dimension(:), allocatable :: nc_ncells_sb ! to store lat / nc_lons real(dp), dimension(:), allocatable :: nc_lons, nc_lats ! time index integer, dimension(:), allocatable :: nc_time_ix ! 2D array for data of one day real(dp), dimension(:, :), allocatable :: nc_var_in #endif contains subroutine input_initialise read(get_config_fid(), INPUT_PARAMETERS) end subroutine input_initialise integer function input_open_file(file_name) result(file_id) ! Abstract input file opener (ascii only) character(len=*), intent(in) :: file_name file_id = open_file(trim(input_dir)//"/"//trim(file_name), 'r') end function input_open_file #ifdef with_netcdf subroutine input_nc_read_climate(flu, humi, mb, ra, subp, tmn, tmx, tx) real(dp), dimension(:), intent(in) :: flu real(dp), dimension(:), intent(out) :: humi integer, intent(in) :: mb real(dp), dimension(:), intent(out) :: ra real(dp), dimension(:), intent(out) :: subp real(dp), dimension(:), intent(out) :: tmn real(dp), dimension(:), intent(out) :: tmx real(dp), dimension(:), intent(out) :: tx ! counters need for day length integer j ! reals need for day length real(dp) :: mm(6) = 0 call input_nc_weighted_mean(1, tx, mb) call input_nc_weighted_mean(2, tmn, mb) call input_nc_weighted_mean(3, tmx, mb) call input_nc_weighted_mean(4, subp, mb) call input_nc_weighted_mean(5, ra, mb) call input_nc_weighted_mean(6, humi, mb) if (NC_DEBUG) then mm = 0 do j = 1, mb mm(1) = mm(1) + tx(j) * flu(j) mm(2) = mm(2) + tmn(j) * flu(j) mm(3) = mm(3) + tmx(j) * flu(j) mm(4) = mm(4) + subp(j) * flu(j) mm(5) = mm(5) + ra(j) * flu(j) mm(6) = mm(6) + humi(j) * flu(j) enddo call log_debug("input_nc_read_climate", & 'ncdf climate (tmean, tmin, tmax, prec, rad, hum) avarage on day:', & int=nc_nday, reals=mm) endif ! update current record number nc_nday = nc_nday + 1 ! If relative humidity is not provided if (isNaN(humi(1)) ) humi = - 999.9 if (isNaN(ra(1)) ) ra = 1. end subroutine input_nc_read_climate subroutine input_nc_weighted_mean(ivar, v_sb, mb) use netcdf, only : nf90_get_var integer, intent(in) :: mb ! variable index integer, intent (in) :: ivar real(dp), intent (out) :: v_sb(:) ! the subbasin variable to be filled integer :: start(3), count(3) ! to supply to nc_get_var ! looping ints integer sb, i ! Read 1 record of NLONS*NLATS*NLVLS values, starting at the beginning ! of the record (the (1, 1, rec) element in the netCDF file). count = (/ nc_nx, nc_ny, 1 /) start = (/ nc_xmin, nc_ymin, nc_nday /) ! read one day/record call input_nc_check_error( nf90_get_var(nc_ids(ivar), nc_var_ids(ivar), nc_var_in, start, count) ) ! apply offset and scale_factor nc_var_in = nc_var_in * nc_attr_vals(ivar, 1) + nc_attr_vals(ivar, 2) ! interpolate and fill v_sb v_sb = 0 do sb = 1, mb do i = 1, nc_ncells_sb(sb) v_sb(sb) = v_sb(sb) + nc_var_in(nc_x_sb(sb, i), nc_y_sb(sb, i)) * nc_weight_sb(sb, i) end do enddo end subroutine input_nc_weighted_mean subroutine input_nc_initialise(iyr, mb) use netcdf, only : nf90_open, nf90_nowrite, nf90_inq_varid integer, intent(in) :: iyr integer, intent(in) :: mb integer :: i logical :: fexists character(len=path_max_length) :: fpaths(6) call log_info('input_nc_initialise', 'Opening netCDF files...') ! allocate arrays with subbasin count allocate(nc_wsum_sb(mb)) nc_wsum_sb = 0 allocate(nc_ncells_sb(mb)) nc_ncells_sb = 0 ! read namelist parameters read(get_config_fid(), nml=NC_CLIMATE_PARAMETERS) ! open nc files, get var ids and compare the dimensions of first file (0), check if others are the same (1) do i = 1, 6 ! check if file exists in input_dir, otherwise interpret as absolute path INQUIRE(file = trim(input_dir) // trim(NC_FNAMES(i)), EXIST = fexists) if (fexists) then fpaths(i) = trim(input_dir) // trim(NC_FNAMES(i)) else fpaths(i) = trim(NC_FNAMES(i)) endif if (NC_DEBUG) call log_debug("input_nc_initialise", trim(fpaths(i))//' ('//trim(NC_VNAMES(i))//')') ! open file, get fileID, check dims and attributes call input_nc_check_error( nf90_open(trim(fpaths(i)), nf90_nowrite, nc_ids(i)) ) call input_nc_check_error( nf90_inq_varid(nc_ids(i), trim(NC_VNAMES(i)), nc_var_ids(i)) ) call input_nc_check_dims(nc_ids(i), i /= 1, trim(NC_VNAMES(i))) call input_nc_check_attr(i) enddo ! check ncgrid.dat, count rows (nc_nrows) and max cell per subbasins (nc_mxc) call input_nc_check_grid(mb) ! allocated arrays with nc_mxc allocate(nc_weight_sb(mb, nc_mxc)) nc_weight_sb = 0 allocate(nc_lon_sb(mb, nc_mxc)) nc_lon_sb = 0 allocate(nc_lat_sb(mb, nc_mxc)) nc_lat_sb = 0 allocate(nc_x_sb(mb, nc_mxc)) nc_x_sb = 0 allocate(nc_y_sb(mb, nc_mxc)) nc_y_sb = 0 ! read ncinfo.dat call input_nc_read_grid(mb) ! convert lon/lat to x/y indecies call input_nc_convert_coordinates(mb) ! find time offset if (NC_REF_YEAR >= 0) then call input_nc_find_time(iyr) ! resets nc_nday else if (NC_OFFSET_DAYS /= 0) then call input_nc_offset_time() ! adds to nc_nday endif ! allocate nc_var_in to read in allocate(nc_var_in(nc_nx, nc_ny)) nc_var_in = 0 ! check for missing values if (NC_DEBUG) call input_nc_check_missing(mb) end subroutine input_nc_initialise subroutine input_nc_check_dims(ncid, icheck, vname) use netcdf, only : & nf90_inq_dimid, nf90_inquire_dimension, nf90_inq_varid, nf90_get_var ! file and variable id to read from integer, intent (in) :: ncid ! check if same as already set values logical, intent (in) :: icheck ! name of variable being checked character(len=*), intent(in) :: vname ! dimension ids to read length integer :: lon_dimid, lat_dimid, rec_dimid ! variable ids to read time, lon/lat values integer :: lon_varid, lat_varid, time_varid integer :: nlo, nla, nre ! Get var ids of lon/lat and time call input_nc_check_error( nf90_inq_dimid(ncid, trim(NC_LON_VNAME), lon_dimid) ) call input_nc_check_error( nf90_inquire_dimension(ncid, lon_dimid, len = nlo) ) call input_nc_check_error( nf90_inq_dimid(ncid, trim(NC_LAT_VNAME), lat_dimid) ) call input_nc_check_error( nf90_inquire_dimension(ncid, lat_dimid, len = nla) ) call input_nc_check_error( nf90_inq_dimid(ncid, trim(NC_TIME_VNAME), rec_dimid) ) call input_nc_check_error( nf90_inquire_dimension(ncid, rec_dimid, len = nre) ) ! report dimensions if (NC_DEBUG) call log_debug("input_nc_check_dims", & 'Number of nc_lons, nc_lats, timesteps:', ints=(/nlo, nla, nre/)) ! check or assign if (icheck) then if (nlo /= nc_nlons) then call log_error("input_nc_check_dims", vname//' has different number of nc_lons', & ints=(/nlo, nc_nlons/)) endif if (nla /= nc_nlats) then call log_error("input_nc_check_dims", vname//' has different number of nc_lats', & ints=(/nla, nc_nlats/)) endif if (nre /= nc_nrecs) then call log_error("input_nc_check_dims", vname//' has different number of time steps', & ints=(/nre, nc_nrecs/)) endif else ! assign nc_nlons = nlo nc_nlats = nla nc_nrecs = nre ! allocate nc_lons, nc_lats, time allocate(nc_lons(nc_nlons)) allocate(nc_lats(nc_nlats)) allocate(nc_time_ix(nc_nrecs)) ! read time, nc_lons and nc_lats call input_nc_check_error( nf90_inq_varid(ncid, trim(NC_LON_VNAME), lon_varid) ) call input_nc_check_error( nf90_get_var(ncid, lon_varid, nc_lons) ) call input_nc_check_error( nf90_inq_varid(ncid, trim(NC_LAT_VNAME), lat_varid) ) call input_nc_check_error( nf90_get_var(ncid, lat_varid, nc_lats) ) call input_nc_check_error( nf90_inq_varid(ncid, trim(NC_TIME_VNAME), time_varid) ) call input_nc_check_error( nf90_get_var(ncid, time_varid, nc_time_ix) ) endif end subroutine input_nc_check_dims subroutine input_nc_check_attr(ivar) use netcdf, only : nf90_get_att, NF90_NOERR ! variable counter integer, intent (in) :: ivar integer :: err, i real(dp) :: attval ! check for variable attributes do i = 1, 4 err = nf90_get_att(nc_ids(ivar), nc_var_ids(ivar), trim(nc_attr_names(i)) , attval) if (err == NF90_NOERR) then nc_attr_vals(ivar, i) = attval ! Attribute default values if not set in addition to initialisation (=0) elseif (i == 1) then ! scale_factor nc_attr_vals(ivar, i) = 1 elseif (i == 4) then ! missing value nc_attr_vals(ivar, i) = nc_attr_vals(ivar, 3) endif if (NC_DEBUG .and. err == NF90_NOERR) then call log_debug("input_nc_check_attr", trim(NC_VNAMES(ivar))//' has a '// & nc_attr_names(i)//' attribute of:', real=nc_attr_vals(ivar, i)) elseif (NC_DEBUG) then call log_debug("input_nc_check_attr", trim(NC_VNAMES(ivar))//' has a '// & nc_attr_names(i)//' default value of:', real=nc_attr_vals(ivar, i)) endif enddo endsubroutine input_nc_check_attr subroutine input_nc_check_grid(mb) integer, intent(in) :: mb integer :: fui, sb integer :: lastsb = 0 integer :: imxc = 0 ! check ncinfo file for length and maxcells fui = input_open_file(NC_GRID) read(fui, *) ! header line do read(fui, *, END=10) sb nc_nrows = nc_nrows + 1 ! count cells per subbasin and update max counter if (sb /= lastsb) imxc = 0 lastsb = sb imxc = imxc + 1 nc_mxc = max(nc_mxc, imxc) ! check subbasin count if (sb > mb) then call log_warn("input_nc_check_grid", "SubbasinID in (line) larger than number of subbasins:", & i1=sb, int=nc_nrows) endif enddo 10 close(fui) ! make sure lastsb==mb if (lastsb /= mb) then call log_warn("input_nc_check_grid", 'Last subbasinID in '//trim(NC_GRID)// & ' does not match the last subbasinID', i1=lastsb, int=mb) endif ! report file length and maxcells per subbasin if (NC_DEBUG) call log_debug("input_nc_check_grid", 'Lines in ncgrid file:', int=nc_nrows) if (NC_DEBUG) call log_debug("input_nc_check_grid", 'Maximum cells per subbasin:', int=nc_mxc) end subroutine input_nc_check_grid subroutine input_nc_read_grid(mb) integer, intent(in) :: mb integer :: fui, i, ii, sb real(dp) :: ilon, ilat, iweight ! open weights and coordinates file fui = input_open_file(NC_GRID) read(fui, *) ! header line do i = 1, nc_nrows ! read line read(fui, *) sb, ilon, ilat, iweight ! count cells nc_ncells_sb(sb) = nc_ncells_sb(sb) + 1 ii = nc_ncells_sb(sb) ! check input and assign call check_range((/ilon/), 'longitude', (/-180., 180./), sb) nc_lon_sb(sb, ii) = ilon call check_range((/ilat/), 'latitude', (/-90., 90./), sb) nc_lat_sb(sb, ii) = ilat call check_range((/iweight/), 'weight', (/0., 1.e6/), sb) nc_weight_sb(sb, ii) = iweight ! weights sum nc_wsum_sb(sb) = nc_wsum_sb(sb) + nc_weight_sb(sb, ii) end do ! close info file again close(fui) ! adjust/make sure weights add up to 1 do sb = 1, mb do ii = 1, nc_ncells_sb(sb) nc_weight_sb(sb, ii) = nc_weight_sb(sb, ii) / nc_wsum_sb(sb) enddo enddo end subroutine input_nc_read_grid subroutine input_nc_convert_coordinates(mb) integer, intent(in) :: mb integer sb, c, i ! loop over subbasins, then cells, then lon/lat stopping at the matching one do sb = 1, mb do c = 1, nc_ncells_sb(sb) do i = 1, nc_nlons if (nc_lons(i) == nc_lon_sb(sb, c)) then nc_x_sb(sb, c) = i nc_xmin = min(nc_xmin, i) nc_xmax = max(nc_xmax, i) exit endif ! if it hasn't found one and got to the end, stop and warn if (i == nc_nlons) then call log_error("input_nc_convert_coordinates", 'Longitude (at subbasin, cell) not in ncdf lon!', & real=nc_lon_sb(sb, c), i1=sb, i2=c) endif end do do i = 1, nc_nlats if (nc_lats(i) == nc_lat_sb(sb, c)) then nc_y_sb(sb, c) = i nc_ymin = min(nc_ymin, i) nc_ymax = max(nc_ymax, i) exit endif ! if it hasn't found one ang got to the end, stop and warn if (i == nc_nlats) then call log_error("input_nc_convert_coordinates", 'Latitude (at subbasin, cell) not in ncdf lat!', & real=nc_lon_sb(sb, c), i1=sb, i2=c) endif end do enddo enddo ! adjust x/y indicies for minimal window nc_nx = nc_xmax - (nc_xmin - 1) nc_ny = nc_ymax - (nc_ymin - 1) do sb = 1, mb do c = 1, nc_ncells_sb(sb) nc_x_sb(sb, c) = nc_x_sb(sb, c) - (nc_xmin - 1) nc_y_sb(sb, c) = nc_y_sb(sb, c) - (nc_ymin - 1) enddo enddo if (NC_DEBUG) call log_debug("input_nc_convert_coordinates", & 'Minimum window size (n nc_lons, n nc_lats):', i1=nc_nx, i2=nc_ny) end subroutine input_nc_convert_coordinates subroutine input_nc_find_time(iyr) integer, intent(in) :: iyr integer ndays_offset, y, d logical :: time_not_found = .true. ! should be replaced by common switch from cod - file integer :: ileap = 1 ndays_offset = (iyr - NC_REF_YEAR) * 365 + NC_OFFSET_DAYS ! add leap years if (ileap == 1) then do y = NC_REF_YEAR, iyr - 1 if ((MOD(y, 4) == 0 .AND. MOD(y, 100) > 0) .OR. MOD(y, 400) == 0) then ndays_offset = ndays_offset + 1 endif enddo endif if (NC_DEBUG) call log_debug("input_nc_find_time", & 'Searching for netCDF time index', int=ndays_offset) do d = 1, nc_nrecs if (ndays_offset == nc_time_ix(d)) then nc_nday = d time_not_found = .false. if (NC_DEBUG) call log_debug("input_nc_find_time", 'Found time index at', int=d) exit ! break from do loop endif enddo ! stop if time index not found if (time_not_found) then call log_error("input_nc_find_time", & 'Cant find time index in climate netCDF file (start, ref. year):', i1=iyr, & i2=NC_REF_YEAR, int=ndays_offset) endif end subroutine input_nc_find_time subroutine input_nc_offset_time() nc_nday = nc_nday + NC_OFFSET_DAYS if (nc_nday <= 0 .OR. nc_nday > nc_nrecs) then call log_error("input_nc_offset_time", & 'Time NC_OFFSET_DAYS leads to out of range index ', int=nc_nday) endif if (NC_DEBUG) call log_error("input_nc_offset_time", & 'Offsetting time index (by) to:', i1=NC_OFFSET_DAYS, int=nc_nday) end subroutine input_nc_offset_time subroutine input_nc_check_missing(mb) use netcdf, only : nf90_get_var integer, intent(in) :: mb integer :: sb, i, ii, nmissing real(dp) :: v call log_debug("input_nc_check_missing", 'Checking for missing values...') nmissing = 0 do i = 1, 6 ! read one day/record call input_nc_check_error( nf90_get_var(nc_ids(i), nc_var_ids(i), nc_var_in, (/ nc_xmin, nc_ymin, nc_nday /), (/ nc_nx, nc_ny, 1 /)) ) ! loop over nc_var_in again and check for missing or fill values do sb = 1, mb do ii = 1, nc_ncells_sb(sb) v = nc_var_in(nc_x_sb(sb, ii), nc_y_sb(sb, ii)) if (v == nc_attr_vals(i, 3) .or. v == nc_attr_vals(i, 4)) then call log_warn("input_nc_check_missing", trim(NC_FNAMES(i))//' ('// & trim(NC_VNAMES(i))//') missing in (subbasin):', i1=sb, & reals=(/nc_lon_sb(sb, ii), nc_lat_sb(sb, ii)/)) nmissing = nmissing + 1 endif enddo enddo enddo if (nmissing > 0) then call log_error("input_nc_check_missing", 'Encountered N missing values:', & int=nmissing) endif endsubroutine input_nc_check_missing subroutine input_nc_close use netcdf, only : nf90_close integer :: i ! deallocate arrays and close open files do i = 1, 6 call input_nc_check_error( nf90_close(nc_ids(i)) ) enddo deallocate(nc_lons) deallocate(nc_lats) deallocate(nc_var_in) endsubroutine input_nc_close subroutine input_nc_check_error(status) use netcdf, only : nf90_noerr, nf90_strerror integer, intent ( in) :: status if (status /= nf90_noerr) then call log_error('input_nc_check_error', trim(nf90_strerror(status))) end if end subroutine input_nc_check_error #endif integer function get_config_fid() integer unit unit = -1 inquire(file=config_file_path, number=unit) if (unit >= 0) then close(unit) end if get_config_fid = open_file(config_file_path) end function get_config_fid function input_count_rows(funit, header) result(nlines) ! Author: stefan.liersch@pik-potsdam.de ! Date: 2009-09-20 ! PURPOSE: counts the number of rows in input file with given file unit ! The file must already be open!!! ! file unit integer, intent(in) :: funit ! true if input file contains a header logical, intent(in) :: header ! counter integer :: nlines ! file status integer :: status ! dummy for reading character(len=read_line_buffer_size) :: a character(len=path_max_length) fname status = 0 nlines = 0 ! rewind file in order to be sure to read the file from the first line rewind(funit) if (header) nlines = nlines - 1 ! count number of rows in file: fname do read(funit, *, IOSTAT=status) a if (status > 0) then ! read error (unlikely since we are only skipping lines) inquire(unit=funit, name=fname) call log_error("input_count_rows", 'Unable to count lines of file '//trim(fname)) else if (status == -1 .or. len_trim(a) == 0) then ! end of file or empty line EXIT end if nlines = nlines + 1 end do rewind(funit) end function input_count_rows subroutine read_logical_column(file_id, column, array, default, index, skip) ! Read column of logical values to fill a 1d array ! ! Optional arguments: ! ------------------- ! column : Name of column to read ! default : Value to fill array if column not found, raises error otherwise ! index : Read column by position/index, since the header is not read ! the column is read from the current line onwards ! skip : Skip rows after header (column parsed) or current position (index) integer, intent(in) :: file_id character(len=*), intent(in), optional :: column logical, dimension(:), intent(out) :: array logical, intent(in), optional :: default integer, intent(in), optional :: index, skip integer i, column_ix character(len=read_line_buffer_size) str_value integer temp, ios if (present(column)) then ! Find column index column_ix = header_column_index(column, file_id) else if (present(index)) then column_ix = index else call log_error("read_logical_column", "column or index argument required.") end if ! Skip lines after the header if (present(skip)) call move_lines(file_id, skip) ! Set default or error if not found if (column_ix == 0 .and. present(default)) then array = default return else if (column_ix == 0) then call input_error_column_not_found(column, file_id) end if ! Fill array do i = 1, size(array) str_value = read_csv_item(column_ix, file_id) read(str_value, *, iostat=ios) temp if (ios /= 0 .or. .not. (temp == 1 .or. temp == 0)) & call input_type_conversion_error(file_id, str_value, "0/1 logical", i+1, & column_ix) array(i) = temp == 1 end do end subroutine read_logical_column subroutine read_real_column(file_id, column, array, default, index, skip, range, closed) ! Read column of real values to fill a 1d array ! ! Optional arguments: ! ------------------- ! column : Name of column to read ! default : Value to fill array if column not found, raises error otherwise ! index : Read column by position/index, since the header is not read ! the column is read from the current line onwards ! skip : Skip rows after header (column parsed) or current position (index) ! range(2) : Check (/lower, upper/) range, column name is required ! closed : Open or closed range, "l"/"u"/"n" lower/upper/neiter, default both integer, intent(in) :: file_id character(len=*), intent(in), optional :: column real(dp), dimension(:), intent(out) :: array real(dp), intent(in), optional :: default, range(2) integer, intent(in), optional :: index, skip character(len=1), intent(in), optional :: closed integer i, column_ix, ios character(len=read_line_buffer_size) str_value character(len=1) clsd if (present(column)) then ! Find column index column_ix = header_column_index(column, file_id) else if (present(index)) then column_ix = index else call log_error("read_real_column", "column or index argument required.") end if ! Skip lines after the header if (present(skip)) call move_lines(file_id, skip) ! Set default or error if not found if (column_ix == 0 .and. present(default)) then array = default return else if (column_ix == 0) then call input_error_column_not_found(column, file_id) end if ! Fill array do i = 1, size(array) str_value = read_csv_item(column_ix, file_id) read(str_value, *, iostat=ios) array(i) if (ios /= 0) & call input_type_conversion_error(file_id, str_value, "real", i+1, column_ix) end do ! Check ranges if (present(range) .and. present(column)) then clsd = "b" if (present(closed)) clsd = closed call check_range(array, trim(column), range, closed=clsd) end if end subroutine read_real_column subroutine read_integer_column(file_id, column, array, default, index, skip, range) ! Read column of integer values to fill a 1d array ! ! Optional arguments: see read_real_column integer, intent(in) :: file_id character(len=*), intent(in), optional :: column integer, dimension(:), intent(out) :: array integer, intent(in), optional :: default, range(2) integer, intent(in), optional :: index, skip integer i, column_ix, ios character(len=read_line_buffer_size) str_value if (present(column)) then ! Find column index column_ix = header_column_index(column, file_id) else if (present(index)) then column_ix = index else call log_error("read_integer_column", "column or index argument required.") end if ! Skip lines after the header if (present(skip)) call move_lines(file_id, skip) ! Set default or error if not found if (column_ix == 0 .and. present(default)) then array = default return else if (column_ix == 0) then call input_error_column_not_found(column, file_id) end if ! Fill array do i = 1, size(array) str_value = read_csv_item(column_ix, file_id) ! Accept fewer items, as newlines at the end of a csv file will break the reader if (len(trim(adjustl(str_value))) .eq. 0) then return end if read(str_value, *, iostat=ios) array(i) if (ios /= 0) & call input_type_conversion_error(file_id, str_value, "integer", i+1, column_ix) end do ! Check ranges if (present(range) .and. present(column)) then call check_int_range(array, trim(column), range) end if end subroutine read_integer_column subroutine read_string_column(file_id, column, array, default, index, skip) ! Read column of strings to fill a 1d array ! ! Optional arguments: see read_real_column integer, intent(in) :: file_id character(len=*), intent(in), optional :: column character(len=*), dimension(:), intent(out) :: array character(len=*), intent(in), optional :: default integer, intent(in), optional :: index, skip integer i, column_ix if (present(column)) then ! Find column index column_ix = header_column_index(column, file_id) else if (present(index)) then column_ix = index else call log_error("read_string_column", "column or index argument required.") end if ! Skip lines after the header if (present(skip)) call move_lines(file_id, skip) ! Set default or error if not found if (column_ix == 0 .and. present(default)) then array = default return else if (column_ix == 0) then call input_error_column_not_found(column, file_id) end if ! Fill array do i = 1, size(array) array(i) = read_csv_item(column_ix, file_id) end do end subroutine read_string_column logical function has_column(column, file_id) character(len=*), intent(in) :: column integer, intent(in) :: file_id has_column = header_column_index(trim(adjustl(column)), file_id) >0 end function has_column integer function header_column_index(column, file_id, required) ! Return index of column in first line of file_id, return 0 if not found character(len=*), intent(in) :: column integer, intent(in) :: file_id logical, intent(in), optional :: required integer iostat, ie character(len=read_line_buffer_size) buffer character(len=read_line_buffer_size*2) scanpart, current_col logical column_not_found rewind(file_id) scanpart = '' column_not_found = .True. header_column_index = 0 iostat = 0 do while (column_not_found .and. iostat == 0) read(file_id, "(a)", advance="no", iostat=iostat) buffer scanpart = trim(adjustl(scanpart))//buffer ie = 1 do while (column_not_found .and. ie > 0) ie = scan(scanpart, ",") if (ie > 0) then current_col = trim(adjustl(scanpart(:ie-1))) scanpart = scanpart(ie+1:) header_column_index = header_column_index + 1 column_not_found = trim(adjustl(current_col)) /= trim(adjustl(column)) else if (iostat < 0) then ! last column current_col = scanpart header_column_index = header_column_index + 1 column_not_found = trim(adjustl(current_col)) /= trim(adjustl(column)) end if end do end do if (column_not_found) then header_column_index = 0 if (present(required)) then if (required) call input_error_column_not_found(column, file_id) end if end if ! Make sure we move to next line if not already if (iostat == 0) call move_lines(file_id, 1) end function header_column_index subroutine move_lines(file_id, n) ! Move n lines forwards or backwards (if negative) in formatted sequential file integer, intent(in) :: file_id, n integer i if (n > 0) then ! move forwards do i = 1, n read(file_id, '()') end do else if (n < 0) then ! move backwards do i = 1, -1*n backspace file_id end do end if end subroutine move_lines function read_csv_item(index, file_id) result(column) ! Return the index'th item from the next line of file_id and move to next line ! Raise error if number of items not found integer, intent(in) :: index, file_id character(len=read_line_buffer_size) :: column character(len=read_line_buffer_size) buffer character(len=read_line_buffer_size*2) scanpart logical column_not_found integer ic, iostat, ie scanpart = '' column_not_found = .True. ic = 0 iostat = 0 do while (column_not_found .and. iostat == 0) read(file_id, "(a)", advance="no", iostat=iostat) buffer scanpart = trim(adjustl(scanpart))//buffer ie = 1 do while (column_not_found .and. ie > 0) ie = scan(scanpart, ",") if (ie > 0) then ! comma found ic = ic + 1 column_not_found = ic /= index if (.not. column_not_found) column = trim(adjustl(scanpart(:ie-1))) scanpart = scanpart(ie+1:) else if (iostat < 0 .and. ic+1 == index) then ! last item column = trim(adjustl(scanpart)) column_not_found = .False. end if end do end do if (column_not_found) then ! TODO error if not found, e.g. because end of file or line else ! Move to next line if not already if (iostat == 0) read(file_id, "(a)") end if end function read_csv_item subroutine input_error_column_not_found(column, file_id) character(len=*), intent(in) :: column integer, intent(in) :: file_id character(len=path_max_length) :: fname inquire(file_id, name=fname) call log_error("input_error_column_not_found", & "Required column "//trim(column)//" not found in "//trim(fname)) end subroutine input_error_column_not_found subroutine input_type_conversion_error(file_id, str, type, line_index, column_index) character(len=*), intent(in) :: str, type integer, intent(in) :: file_id, line_index, column_index character(len=path_max_length) :: fname inquire(file_id, name=fname) call log_error("input_type_conversion_error", "Failed to convert "//trim(str)// & " to "//trim(type)//" of "//trim(fname)//" in (line, column)", i1=line_index, i2=column_index) end subroutine input_type_conversion_error subroutine print_splash(version) use utilities, only: colourise character(len=*), intent(in) :: version integer :: nlsp, nrsp character(len=12) :: l1 character(len=14) :: l2, l4 character(len=16) :: l3 l1 = colourise("~~~", 'yellow') l2 = colourise("~~~~~", 'green') l3 = colourise("~~~~~~~", 'cyan') l4 = colourise("~~~~~", 'blue') nlsp = floor((49 - len(version)) / 2.) nrsp = ceiling((49 - len(version)) / 2.) write(*, "(A)") "", & " + . . . . . . . . . . . . . . . . . . . . . . . . . . +", & " . .", & " . . //////// \\ \\ // || || || .", & " . / \ \\ \\ \\ // || ||\\//|| .", & " . /"//l1//"\ \\\\\\\ \\ //\\ // || || \/ || .", & " . /"//l2//"\ // \\// \\// || || || .", & " . (" //l3// ") /////// \\ \\ || || || .", & " . \"//l4//"/ -------------------------------------- .", & " . ----- Soil and Water Integrated Model .", & " . .", & " ."//repeat(" ", nlsp)//"* "//version//" *"//repeat(" ", nrsp)//".", & " + . . . . . . . . . . . . . . . . . . . . . . . . . . +", "" end subroutine print_splash subroutine print_help(errormsg) use utilities, only: log_error character(len=*), optional :: errormsg write(*, "(A)") & "usage: swim [options] [parameter-nml]", & "", & " parameter-nml path to the parameter namelist", & "", & "options:", & " -h, --help show this help message", & " -v, --version print version", & " -d, --defaults [module] print default parameters", & " -o, --output-variables print output variables (requires parameter-nml)" ! End program either with error or silently if (present(errormsg)) then write(*, *) ! empty line btw help and error call log_error("print_help", errormsg) end if stop end subroutine print_help subroutine parse_commandline_arguments(version, print_output_variables, print_defaults) character(len=*), intent(in) :: version logical, intent(out) :: print_output_variables character(len=*), intent(out) :: print_defaults character(len=path_max_length) :: arg integer i, iarg ! Number of commandline arguments iarg = iargc() ! without arguments only show help if (iarg == 0) call print_help("parameter-nml is missing") print_output_variables = .False. print_defaults = " " do i = 1, iarg call getarg(i, arg) ! Options if (arg(:1) == "-") then ! remove first - if long form if (arg(1:2) == "--") arg = arg(2:) select case (trim(arg(2:))) case ("h", "help") call print_splash(version) call print_help case ("d", "defaults") ! Interpret next argument as module call getarg(i + 1, print_defaults) if (trim(print_defaults) == "") print_defaults = "all" case ("o", "output-variables") print_output_variables = .True. case ("v", "version") write(*, "(A)") version stop case default call print_help("Unknown option: "//trim(arg)) end select else ! Positional argument config_file_path = arg end if end do end subroutine parse_commandline_arguments end module input