! !     FILE utilities.f95
!

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! INDICES:
!      i, ida - day
!      j - subbasin
!      jea, je - hydrotope or HRU
!      k - soil type
!      n - land use type
!      l - soil layers
!
! PHASE VARIABLES & FLUXES:
! HYDROLOGY:
!      swe(mb, meap) = soil water content, mm
!      ste(mb, meap, ml) = water storage in a layer, mm
!      snoa(mb, meap) = snow water content, mm
!      flate(mb, meap, ml) = lateral subsurface flow for layer, mm
!      poe(mb, meap, ml) = percolation for layer, mm
! CROP:
!      cva(mb, meap) = vegetation cover, kg/ha
!      dm(mb, meap) = total biomass, kg/ha
!      rsd(mb, meap) = crop residue, kg/ha
! NUTRIENTS:
!      ano3(mb, meap, ml) = nitrate (NO3-N) content in a layer, kg/ha
!      anora(mb, meap, ml) = active org. N content in a layer, kg/ha
!      anors(mb, meap, ml) = stable org. N content in a layer, kg/ha
!      plab(mb, meap, ml) = labile P content in a layer, kg/ha
!      porg(mb, meap, ml) = organic P in a layer, kg/ha
!      pma(mb, meap, ml) = active mineral P in a layer, kg/ha
!      pms(mb, meap, ml) = stable mineral P in a layer, kg/ha
!      fop(mb, meap) = fresh org N, kg/ha
!      fon(mb, meap) = fresh org P, kg/ha
! OTHER IMPORTANT VARIABLES:
!      mstruc(mb, meap, 4) = HRU structure vector. land use, soil, wet, manag
!      frar(mb, meap) = fractional area of hydrotope in subbasin
!      te(mb, meap, ml) = daily ave temp at the bottom of layer, degree !
!      te0(mb, meap) = bare soil surface temp, degree !
!      smx(mb, meap) = retention factor, corresponding cn1
!      wf(2, mb, meap) = shape parameters for calc. of retention
!      igro(mb, meap) = vegetation index, =1 if vegetation is growing
!      g(mb, meap) = fraction of heat units to maturity accumulated
!      alai(mb, meap) = leaf area index
!      rwt(mb, meap) = fraction of root weight
!      hia(mb, meap) = harvest index
!      hiad(mb, meap) = harvest index, adjusted
!      canstor(mb, meap) = canopy water storage, mm
!      ws(mb, meap) = water stress factor
!      huharv(mb, meap) = harvest index heat unit
!      rd(mb, meap) = root depth, mm
!      yld(mb, meap) = crop yield, kg/ha
!      cklsp(mb, meap) = combined c, k, ls, p factor
!      snup(mb, meap) = N uptake, kg/ha
!      spup(mb, meap) = P uptake, kg/ha
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
module utilities

  implicit none

  ! Global real type (inspired by F90 netCDF)
  integer, parameter :: dp = 4 !selected_real_kind(P = 13, R = 307)

    ! i/o constants
  integer, parameter :: &
    identifier_max_length = 32, &
    ! Should be the number of registered variables, this must be a parameter
    output_max_variables = 32, &
    path_max_length = 1024, &
    ! Max length of output labels
    output_label_max_length = 64

  ! These variables should renamin in utilities
  real(dp), save :: abfzwo
  real(dp), save :: del0
  real(dp), save :: evafac
  real(dp), save :: gwafac

  integer :: ifod = 59
  integer :: ifoe = 60
  integer :: ifom = 61
  integer :: ihei = 50
  integer :: imai = 25
  integer :: imea = 46
  integer :: ipas = 47
  integer :: irp = 40
  integer :: iwet = 55
  integer :: iwetf = 63
  integer, save, dimension(4) :: k10 = (/ 73, 24, 881, 52/)
  integer, save, dimension(4) :: k11 = (/ 6, 302, 597, 3/)
  real(dp), save :: ndmo(12) !
  character(len = 4) :: prog
  ! rotated crop (po, ma, sb, wr, wb, wwII, wrape, ww) see initcrp_xx.f
  character(len = 5) :: rot_variant
  character(len = 8) :: till
  ! State for gammad utility
  real(dp), save :: xn1 = 10.

  ! Log level only useful for debugging
  integer, parameter :: log_debug_level   = 10
  ! Log level generated by the normal execution of the program.
  integer, parameter :: log_info_level    = 20
  ! Log level which indicates parts of the program arnt behaving as it ideally should.
  integer, parameter :: log_warning_level = 30
  ! Log level which indicates an error has occurred fatally ending execution.
  integer, parameter :: log_error_level   = 40
  ! progress reminder
  integer :: log_in_progess = 0

  type :: logger
    ! An object to handle output of information about the executing
    ! program to the terminal and to a log-file.

    ! Unit corresponding to STDOUT
    integer                       :: stdout = 6
    ! Unit corresponding to STDERR
    integer                       :: stderr = 0
    ! Unit corresponding to log-file
    integer                       :: fileunit = -1
    ! Name of the log-file
    character(len=path_max_length) :: logfile
    ! Cutoff for which messages will be written to STDERR.
    integer                       :: stderr_threshold = log_debug_level
    ! Cutoff for which messages will be written to STDOUT.
    integer                       :: stdout_threshold = log_debug_level
    ! Cutoff for which messages will be written to the log-file.
    integer                       :: logfile_threshold = huge(1)
  end type logger

  ! The main logger object for a program to use
  type(logger), save :: master_log

contains

  integer function days_in_month(month, year)

    integer, intent(in) :: month, year
    integer, dimension(12) :: ndays, ndays_leap

    days_in_month = 0
    ndays = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
    ndays_leap = (/ 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
    if (is_leap_year(year) ) then
      days_in_month = ndays_leap(month)
    else
      days_in_month = ndays(month)
    end if

  end function days_in_month

  integer function day_of_month(ida, month, year, nc)

    integer, intent(in) :: ida
    integer, dimension(13), intent(in) :: nc
    integer, intent(in) :: month, year
    integer :: n

    ! ida = current day (day of year) used in SWIM as global variable
    ! nc(month) = number of days passed in the beginning of month (common.f90) global variable in SWIM
    n = nc(month)
    if (.NOT.is_leap_year(year) .AND. month > 2 ) n = n - 1
    day_of_month = ida - n

  end function day_of_month

  logical function is_leap_year(year)

    ! returns .true. if year is a leap year, otherwise .false.
    integer, intent(in) :: year

    if ((MOD(year, 4) == 0 .AND. MOD(year, 100) > 0) .OR. MOD(year, 400) == 0) then
      is_leap_year = .true.
    else
      is_leap_year = .false.
    end if

  end function is_leap_year

  function hydrograph_storage_location(icodes, ihouts, inum1s, mhyd, nsubs, subnr)
    ! Function returns an array with the hydrograph storage locations (.fig file)
    ! for subbasins listed in subnr

    ! number of subbasins asking for ihout
    integer, intent(in) :: nsubs
    ! 1 - dim array containing subbasin numbers
    integer , intent(in) :: subnr(nsubs)
    integer, dimension(:), intent(in) :: icodes
    integer, dimension(:), intent(in) :: ihouts
    integer, dimension(:), intent(in) :: inum1s
    integer, intent(in) :: mhyd
    ! array of hydrograph storage locations
    integer :: hydrograph_storage_location(nsubs)
    integer :: i, j

    hydrograph_storage_location = 0

    ! identify hydrograph storage location for each output gauge (subbasin)
    do i = 1, mhyd
      if (icodes(i) == 2 ) then ! ROUTE command
        do j = 1, nsubs
          if (inum1s(i) == subnr(j) ) then
            ! The desired hydrograph storage location is always the next (+1)
            !ihouts_out(j) = ihouts(i) + 1
            hydrograph_storage_location(j) = ihouts(i) + 1
          end if
        end do
      end if
    end do

    ! Some subbasins will certainly not be captured by the algorithm above.
    ! Normally these are headwater subbasins.
    ! The following loop captures those headwater subbasins.
    do j = 1, nsubs
      if (hydrograph_storage_location(j) == 0) then
        hydrograph_storage_location(j) = subnr(j)
      end if
    end do

  end function hydrograph_storage_location

  function gamma_distribution(rn1, ai, k7)
    !**** PURPOSE: THIS FUNCTION PROVIDES NUMBERS rn1 FROM A GAMMA DISTRIBUTION,
    !              GIVEN TWO RANDOM NUMBERS ai & k7
    !**** CALLED IN: ALPHA
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> STATIC PARAMETERS
    !      fu = local par
    !      rn = local par
    !      xn1 = local par
    !      xx = local par
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    integer k7
    real(dp) rn1, ai, fu, rn, xx, gamma_distribution

    dimension k7(4)

    do while (.true.)
      gamma_distribution = rn1
      xx = rn1 * ai
      rn = random_n(k7)
      fu = xx ** xn1 * exp(xn1 * (1. - xx))
      rn1 = rn
      if (fu .ge. rn) exit
    end do
    return
  end function gamma_distribution

  function random_n(k)
    !**** PURPOSE: THIS FUNCTION PROVIDES RANDOM NUMBERS RANGING FROM 0. TO 1.
    !**** CALLED IN: GAMMAD, INIT, READSUB
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> STATIC
    !      i = local par
    !      >>>>>
    !calc randn
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    integer i, k
    real(dp) random_n

    dimension k(4)
    k(4) = 3 * k(4) + k(2)
    k(3) = 3 * k(3) + k(1)
    k(2) = 3 * k(2)
    k(1) = 3 * k(1)
    i = k(1) / 1000
    k(1) = k(1) - i * 1000
    k(2) = k(2) + i
    i = k(2) / 100
    k(2) = k(2) - 100 * i
    k(3) = k(3) + i
    i = k(3) / 1000
    k(3) = k(3) - i * 1000
    k(4) = k(4) + i
    i = k(4) / 100
    k(4) = k(4) - 100 * i
    random_n = (((float(k(1)) * .001 + float(k(2))) * .01 + float(k(3))) * .001 + &
        float(k(4))) * .01
    return
  end function random_n

  function indeces_in_string_array(entries, array)  result(outarray)
    ! Get the indeces in array of entries
    ! outarray must be the same length as entries
    character(len=*), dimension(:), intent(in) :: entries, array
    integer, dimension(size(entries)) :: outarray
    integer i

    do i = 1, size(entries)
      outarray(i) = string_index(entries(i), array)
    end do
  end function indeces_in_string_array

  function int_index(id, array)  result(idx)
    ! Get the indeces in array of entries, array must only contain unique
    ! values otherwise only the indeces of the first occurence is returned
    ! outarray must be the same length as entries
    integer, intent(in) :: id
    integer, dimension(:), intent(in) :: array
    integer :: idx
    integer i

    idx = 0
    do i = 1, size(array)
      if (array(i) == id) then
        idx = i
        exit
      end if
    end do
    if (idx == 0) then
      call log_error("indeces_in_int_array", "ID not in array:", int=id, ints=array)
    end if
  end function int_index

  function indeces_in_int_array(entries, array)  result(outarray)
    ! Get the indeces in array of entries, array must only contain unique
    ! values otherwise only the indeces of the first occurence is returned
    ! outarray must be the same length as entries
    integer, dimension(:), intent(in) :: entries, array
    integer, dimension(size(entries)) :: outarray
    integer i

    do i = 1, size(entries)
      outarray(i) = int_index(entries(i), array)
    end do
  end function indeces_in_int_array

  function date_time_str()
    ! Returns the formatted current date and time.

    character(len=19) :: date_time_str
    integer, dimension(8) :: time_vals
    character(len=44), parameter :: time_format = '(i4,"-",i2,"-",i2,1x,'// &
                                                  'i2.2,":",i2.2,":",i2.2)'
    call date_and_time(values=time_vals)
    write(date_time_str, time_format) time_vals(1:3), time_vals(5:7)
  end function date_time_str

  function colourise(str, colour) result(cstr)
    ! Add ANSI(8) colour codes around the string.
    character(len=*), intent(in) :: str, colour
    character(len=len(str)+9) :: cstr
    character(len=1), parameter :: c_esc = achar(27)
    character(len=2), parameter :: c_start = c_esc // '['
    character(len=1), parameter :: c_end = 'm'
    character(len=*), parameter :: c_clear = c_start // '0' // c_end
    character(len=*), parameter, dimension(2, 6) :: codes = reshape((/ &
       'red    ', '31     ', &
       'green  ', '32     ', &
       'yellow ', '33     ', &
       'blue   ', '34     ', &
       'magenta', '35     ', &
       'cyan   ', '36     ' &
       /), (/2, 6/))

    cstr = c_start // trim(codes(2, string_index(colour, codes(1, :)))) // c_end // &
      trim(adjustl(str)) // c_clear
  end function colourise

  integer function log_str2level(str)
    ! Convert a str log level to a level number
    character(len=*), intent(in) :: str
    select case (str)
        case ('d', 'debug')
            log_str2level = log_debug_level
        case ('i', 'info')
            log_str2level = log_info_level
        case ('w', 'warn', 'warning')
            log_str2level = log_warning_level
        case ('e', 'error')
            log_str2level = log_error_level
        case default
            log_str2level = log_debug_level
    end select
  end function log_str2level

  function log_create(logfile, stderr_level, stdout_level, logfile_level) result(log)
    type(logger) :: log
    ! Name of the log-file to which output will be written
    character(len=*), intent(in)  :: logfile
    ! Threshold priority, at and above which messages will be written to standard error.
    character(len=*), intent(in), optional :: stderr_level
    ! Threshold priority, at and above which messages will be written to standard out.
    character(len=*), intent(in), optional :: stdout_level
    ! Threshold priority, at and above which messages will be written to the log file.
    character(len=*), intent(in), optional :: logfile_level

    log%logfile = logfile
    if (present(stderr_level)) log%stderr_threshold = log_str2level(stderr_level)
    if (present(stdout_level)) log%stdout_threshold = log_str2level(stdout_level)
    if (present(logfile_level)) log%logfile_threshold = log_str2level(logfile_level)
    ! open file if needed
    if (log%logfile_threshold <= log_error_level) log%fileunit = open_file(logfile, 'w')
  end function log_create

  function log_format_message(level, tag, message, colour, datetime) result(msg)
    ! Format log message with tag, message and optionally coloured level

    ! The level for which to get the designator
    integer, intent(in) :: level
    character(len=*), intent(in) :: tag
    character(len=*), intent(in) :: message
    ! Whether to colourise the level tag, should only be done if the designator is
    ! printed to the terminal. Defaults to `.true.`.
    ! Whether to add datetime stamp at beginning, defaults to .true.
    logical, intent(in), optional :: colour, datetime
    character(len=23+len(tag)+2+21+len(message)) :: msg
    character(len=21) :: levelstr, dt
    character(len=len(tag)+2) :: tagb
    integer :: descriminator
    logical :: col

    if (present(colour)) then
      col = colour
    else
      col = .true.
    end if
    descriminator = level/10
    if (descriminator < 1) then
        levelstr = ''
    else if (descriminator < 2) then
        levelstr = '<DEBUG>'
      if (col) levelstr = colourise(levelstr, 'cyan')
    else if (descriminator < 3) then
        levelstr = '<INFO>'
      if (col) levelstr = colourise(levelstr, 'blue')
    else if (descriminator < 4) then
        levelstr = '<WARNING>'
      if (col) levelstr = colourise(levelstr, 'yellow')
    else if (descriminator < 5) then
        levelstr = '<ERROR>'
      if (col) levelstr = colourise(levelstr, 'red')
    end if
    dt = ''
    if (present(datetime)) then
      if (datetime) dt = date_time_str()
    else
      dt = "["//date_time_str()//"]"
    end if
    tagb = tag
    if (len_trim(tag) > 0) tagb = '['//tag//']'
    ! put together
    write(msg, "(3a,(1x,a))") trim(levelstr), trim(dt), trim(tagb), trim(message)
  end function log_format_message

  subroutine log_write(fid, msg, i1, i2, int, real, ints, reals, advance, clear)
    ! Write msg to fid appending various numerical values
    integer, intent(in) :: fid
    character(len=*), intent(in) :: msg
    ! optional (indeces), single int/real or 1D arrays to append to message
    integer, intent(in), optional :: i1, i2, int, ints(:)
    real(dp), intent(in), optional :: real, reals(:)
    ! add line break (default true)
    logical, intent(in), optional :: advance
    ! clear n characters before writing (default 0)
    integer, intent(in), optional :: clear

    integer i

    ! clear line if needed
    if (present(clear)) then
      if (clear > 0) write(fid, '(60a)', advance='no') (char(8), i=1, clear)
    end if

    write(fid, '(a)', advance='no') trim(msg)
    if (present(i1) .and. .not. present(i2)) &
      write(fid, '(1x,"(",a,")")', advance='no') trim(to_string(int=i1))
    if (present(i2) .and. .not. present(i1)) &
      write(fid, '(1x,"(",a,")")', advance='no') trim(to_string(int=i2))
    if (present(i2) .and. present(i1)) &
      write(fid, '(1x,"(",a,", ",a,")")', advance='no') &
      trim(to_string(int=i1)), trim(to_string(int=i2))
    if (present(int)) write(fid, '(1x,a)', advance='no') trim(to_string(int=int))
    if (present(real)) write(fid, '(1x,a)', advance='no') trim(to_string(real=real))
    if (present(ints)) write(fid, '(*(1x,a))', advance='no') &
      (trim(to_string(int=ints(i))), i=1, size(ints))
    if (present(reals)) write(fid, '(*(1x,a))', advance='no') &
      (trim(to_string(real=reals(i))), i=1, size(reals))
    ! New line unless advance=False
    if (present(advance)) then
      if (advance) write(fid, '(a)')
    else
      write(fid, '(a)')
    end if
  end subroutine log_write

  subroutine log_message(source, level, message, log, i1, i2, int, real, ints, reals)
    ! Write the provided message to STDERR, STDOUT, and/or a log-file, based on level.

    ! The name of the procedure which produced the error
    character(len=*), intent(in) :: source
    ! The importance of the message, determining where it will be written.
    integer, intent(in)          :: level
    character(len=*), intent(in) :: message
    ! optional logger, defaults to master_log
    type(logger), intent(in), optional :: log
    ! optional index, single int/real or 1D arrays to append to message
    integer, intent(in), optional :: i1, i2, int, ints(:)
    real(dp), intent(in), optional :: real, reals(:)

    type(logger) :: lg
    integer :: cl = 0

    if (present(log)) then
        lg = log
    else
        lg = master_log
    end if
    if (level >= lg%stderr_threshold) then
      if (lg%stderr == lg%stdout) cl = log_in_progess
      call log_write(lg%stderr, &
        log_format_message(level, source, message, datetime=.false.), &
        i1, i2, int, real, ints, reals, clear=cl)
    else if (level >= lg%stdout_threshold) then
      call log_write(lg%stdout, &
        log_format_message(level, source, message, datetime=.false.), &
        i1, i2, int, real, ints, reals, clear=log_in_progess)
    end if
    if (lg%fileunit > 0 .and. level >= lg%logfile_threshold) then
      call log_write(lg%fileunit, &
        log_format_message(level, source, message, colour=.false.), &
        i1, i2, int, real, ints, reals)
    end if
  end subroutine log_message

  subroutine log_progress(source, step, total)
    ! Updates progress bar on stdout as INFO
    character(len=*), intent(in) :: source
    integer, intent(in) :: step, total
    character(len=58), parameter :: &
      bar = '==========================================================', &
      emp = '                                                          '
    integer clen, length

    ! Only print progress bar if stdout == info
    if (master_log%stdout_threshold /= log_info_level) return

    ! Length of bar - 2*[] - <INFO>
    length = len(bar) - len_trim(source) - 2 - 7 - 2
    clen = int(length*step/total)
    call log_write(master_log%stdout, &
      log_format_message(20, source, "["//bar(:clen)//emp(len(emp)-length+clen+1:)//"]", datetime=.false.), &
      advance=(step == total), clear=len(bar))
    ! remember that in progress
    if (step /= total) then
      log_in_progess = len(bar)
    else
      log_in_progess = 0
    end if
  end subroutine log_progress

  subroutine log_debug(source, message, log, i1, i2, int, real, ints, reals)
    ! Writes log_debug_level information to STDERR, STDOUT, and/or a log-file.

    ! The name of the procedure which produced the error
    character(len=*), intent(in) :: source
    ! The information to be written.
    character(len=*), intent(in) :: message
    ! optional logger, defaults to master_log
    type(logger), intent(in), optional :: log
    ! optional index, single int/real or 1D arrays to append to message
    integer, intent(in), optional :: i1, i2, int, ints(:)
    real(dp), intent(in), optional :: real, reals(:)
    call log_message(source, log_debug_level, message, log, i1, i2, int, real, ints, reals)
  end subroutine log_debug

  subroutine log_info(source, message, log, i1, i2, int, real, ints, reals)
    ! Writes run-time information to STDERR, STDOUT, and/or a log-file.

    ! The name of the procedure which produced the error
    character(len=*), intent(in) :: source
    ! The information to be written.
    character(len=*), intent(in) :: message
    ! optional logger, defaults to master_log
    type(logger), intent(in), optional :: log
    ! optional index, single int/real or 1D arrays to append to message
    integer, intent(in), optional :: i1, i2, int, ints(:)
    real(dp), intent(in), optional :: real, reals(:)
    call log_message(source, log_info_level, message, log, i1, i2, int, real, ints, reals)
  end subroutine log_info

  subroutine log_warn(source, message, log, i1, i2, int, real, ints, reals)
    ! Writes warning information to STDERR, STDOUT, and/or a log-file.

    ! The name of the procedure which produced the error
    character(len=*), intent(in) :: source
    ! The information to be written.
    character(len=*), intent(in) :: message
    ! optional logger, defaults to master_log
    type(logger), intent(in), optional :: log
    ! optional index, single int/real or 1D arrays to append to message
    integer, intent(in), optional :: i1, i2, int, ints(:)
    real(dp), intent(in), optional :: real, reals(:)
    call log_message(source, log_warning_level, message, log, i1, i2, int, real, ints, reals)
  end subroutine log_warn

  subroutine log_error(source, message, log, i1, i2, int, real, ints, reals)
    ! Writes error information to STDERR, STDOUT, and/or a log-file.

    ! The name of the procedure which produced the error
    character(len=*), intent(in) :: source
    ! The information to be written.
    character(len=*), intent(in) :: message
    ! optional logger, defaults to master_log
    type(logger), intent(in), optional :: log
    ! optional index, single int/real or 1D arrays to append to message
    integer, intent(in), optional :: i1, i2, int, ints(:)
    real(dp), intent(in), optional :: real, reals(:)
    call log_message(source, log_error_level, message, log, i1, i2, int, real, ints, reals)
    ! close log and stop execution
    if (log%fileunit > 0) close(log%fileunit)
    call abort
  end subroutine log_error

  integer function open_file(fname, status)
    ! This function opens a file and returns the unit id (a free unit number between 21-1000,
    ! if file open was successful. The file status remains open!
    ! If file opening was not successful, the functions stops the programme and closes the file.
    character(len=*), intent(in) :: fname
    ! r for reading, w for writing
    character(len=1), intent(in), optional :: status
    character(len=1) :: st
    integer :: opensucc
    integer ui, iostat
    logical opened

    st = 'r'
    if (present(status)) st = status
    ! get free unit id
    do ui = 21, 1000
      inquire(unit = ui, opened = opened, iostat = iostat)
      if (iostat .ne. 0) cycle
      if (.not. opened) exit
    end do
    open_file = ui
    ! open file
    if (trim(st) == 'r') then
      open(ui, file=trim(fname), status='OLD', IOSTAT=opensucc)
    else if (trim(st) == "w") then
      open(ui, file=trim(fname), IOSTAT=opensucc)
    endif
    if (opensucc /= 0) then
      close(ui)
      call log_error("open_file", "ERROR while opening file: "//trim(fname)//" does it exist?")
    end if
  end function open_file

  subroutine extend_unique_string(nonunique, unique)
    ! Add entries from nonunique to the next blank of unique unless already present
    character(len=*), dimension(:), intent(in) :: nonunique
    character(len=*), dimension(:), intent(inout) :: unique
    logical add
    integer i, ii, iv
    iv = 1
    ! Shift iv counter to allow for already present unique entries
    do ii = 1, size(unique)
      if (trim(unique(ii)) == "") then
        exit
      end if
      iv = iv + 1
    end do
    ! Add nonunique to unique unless already in unique
    do i = 1, size(nonunique)
      add = .True.
      do ii = 1, size(unique)
        if (trim(nonunique(i)) == trim(unique(ii))) then
          add = .False.
          exit
        else if (trim(unique(ii)) == "") then
          exit
        end if
      end do
      if (add) then
        if (iv > size(unique)) then
          call log_error('extend_unique_string', &
            "Trying to extend unique array autside size.", i1=iv, int=size(unique))
        end if
        unique(iv) = nonunique(i)
        iv = iv + 1
      end if
    end do
  end subroutine extend_unique_string

  integer function string_index(str, array)
  ! Return index of first occurence of str in array, return 0 if not in
  character(len=*), intent(in) :: str, array(:)
  integer i
  string_index = 0
  do i = 1, size(array)
    if (trim(array(i)) == trim(str)) then
      string_index = i
      exit
    end if
  end do
  end function string_index

  subroutine input_split(str, delim, before)
    !     Routine splits the string 'str' by the first occurence of delim.
    !     The characters before the found delimiter are output in 'before'.
    !     The characters after the found delimiter are output in 'str'.


    character(len = 300) str
    character(len = 200) before
    character delim
    integer ipos

    str = trim(str)
    if (len(str) == 0) return ! string str is empty

    ipos = index(str, delim)
    if (ipos > 0) then ! next character is a delimiter
      before = str(1:ipos - 1)
      str = str(ipos + 1:)
      before = trim(before)
      str = trim(str)
    end if
    return

  end subroutine input_split

  character(len=64) function to_string(real, int, logical, fmt)
    ! Convert any real/int/logical value to a left-aligned string using default
    ! formatting. Since functions cant return variable length strings, this
    ! should likely be combined with trim().
    ! Examples:
    !  trim(to_string(1.33))
    !  trim(to_sting(int=1200))
    real(dp), intent(in), optional :: real
    integer, intent(in), optional :: int
    logical, intent(in), optional :: logical
    character(len=*), intent(in), optional :: fmt

    character(len=len(to_string)) :: number, decim, rfmt
    integer i

    to_string = ""
    number = ""
    if (present(real)) then
      if (present(fmt)) then
        rfmt = "("//fmt//")"
      else
        ! Find last significant decimal place
        do i=0, 16
          if (mod(real*10**i, 1.) == 0) exit
        end do
        write(decim, *) i
        rfmt = '(F64.'//trim(adjustl(decim))//')'
      end if
      write(number, rfmt) real
    else if (present(int)) then
      if (present(fmt)) then
        write(number, "("//fmt//")") int
      else
        write(number, *) int
      end if
    else if (present(logical)) then
      write(number, *) logical
    end if
    to_string = adjustl(number)
  end function to_string

  subroutine check_range(values, name, range, index, closed)
    ! Check if the real array values is within (/lower, upper/) range and terminate
    ! with name if not. values must be a 1D array but a single
    ! value can also be checked by passing (/value/) and an optional index to report.
    ! The closed optional argument specifies the inclusion/exclusion of the bounds:
    ! closed="l"/"u"/"n" means lower/upper/neither is included, default is both.
    ! Examples:
    !  call check_range((/1., 2., 3./), "array", (/0., 2./))
    !  call check_range((/2./), "val", (/0., 2./), closed="n")
    !  call check_range((/3./), "array", (/0., 2./), index=100)
    real(dp), intent(in) :: values(:)
    character(*), intent(in) :: name
    real(dp), intent(in) :: range(2)
    integer, intent(in), optional :: index
    character(len=1), intent(in), optional :: closed

    integer i, nv, counter
    real(dp) va
    logical clsd(2)

    if (range(1) >= range(2)) &
      call log_error("check_range", "range is not (/lower, upper/) values", reals=range)
    clsd = .True.
    if (present(closed)) then
      if (closed == "n") clsd = .False.
      if (closed == "l") clsd(2) = .False.
      if (closed == "u") clsd(1) = .False.
    end if
    counter = 0
    nv = size(values)
    do i = 1, nv
      va = values(i)
      ! check range depending on range closure
      if ((clsd(1) .and. va < range(1)) .or. (.not.clsd(1) .and. va <= range(1)) .or. &
          (clsd(2) .and. va > range(2)) .or. (.not.clsd(2) .and. va >= range(2))) then
          if (nv > 1) &
            call log_warn('check_range', name, i1=i, real=va)
          counter = counter + 1
      endif
    end do
    ! Array with counted out of range values, indeces reported above
    if (counter > 0 .and. nv > 1) call out_of_range_error(name, range, counter=counter)
    ! Single value with index
    if (counter == 1 .and. present(index)) &
      call out_of_range_error(name, range, va, index=index)
    ! Single value without index
    if (counter == 1) call out_of_range_error(name, range, values(1))
  end subroutine check_range

  subroutine check_int_range(values, name, range, index)
    ! Check if the integer array values is within (/lower, upper/) range and terminate
    ! with name if not. values must be a 1D array but a single
    ! value can also be checked by passing (/value/) and an optional index to report.
    ! Examples:
    !  call check_int_range((/1, 2, 3/), "array", (/0, 2/))
    !  call check_int_range((/3/), "val", (/0, 2/))
    !  call check_int_range((/3/), "array", (/0, 2/), index=100)
    integer, intent(in) :: values(:)
    character(*), intent(in) :: name
    integer, intent(in) :: range(2)
    integer, intent(in), optional :: index

    integer i, nv, counter

    if (range(1) >= range(2)) &
      call log_error("check_int_range", "range is not (/lower, upper/) values", &
        ints=range)
    counter = 0
    nv = size(values)
    do i = 1, nv
      if (values(i) < range(1) .or. values(i) > range(2)) then
          ! Only print value if not single cell array
          if (nv > 1) call log_warn('check_int_range', name, i1=i, int=values(i))
          counter = counter + 1
      endif
    end do
    ! Array with counted out of range values, indeces reported above
    if (counter > 0 .and. nv > 1) call out_of_range_error(name, int_range=range, counter=counter)
    ! Single value with index
    if (counter == 1 .and. present(index)) &
      call out_of_range_error(name, int_range=range, int_value=values(1), index=index)
    ! Single value without index
    if (counter == 1) call out_of_range_error(name, int_range=range, int_value=values(1))
  end subroutine check_int_range

  subroutine out_of_range_error(name, range, value, int_range, int_value, counter, index)
    ! Abstracted fatal error used in check_range and check_int_range
    character(*), intent(in) :: name
    real(dp), intent(in), optional :: value, range(2)
    integer, intent(in), optional :: int_value, int_range(2), counter, index

    character(len=64) strvalue
    character(len=2*len(strvalue)) strrange, strat

    strrange = ""
    strvalue = ""
    if (present(range)) then
      if (present(value)) strvalue = to_string(value)
      strrange = trim(to_string(range(1)))//" - "//to_string(range(2))
    else if (present(int_range)) then
      if (present(int_value)) strvalue = to_string(int=int_value)
      strrange = trim(to_string(int=int_range(1)))//" - "//to_string(int=int_range(2))
    else
      call log_error("out_of_range_error", "range or int_range must be given")
    end if
    ! Array with counter
    if (present(counter)) then
      call log_error("out_of_range_error", "Number of values of "//name//" are out "// &
        "of range ("//trim(strrange)//")", int=counter)
    ! Single value with/without index
    else
      strat = ""
      if (present(index)) strat = "("//trim(to_string(int=index))//")"
      call log_error("out_of_range_error", name//trim(strat)//" = "//trim(strvalue)// &
        " is out of range", reals=range)
    end if
  end subroutine out_of_range_error

end module utilities