Fortran plot

module plotutil
  implicit none

  integer, parameter :: rk = 8

  character(3000) :: line

contains
  subroutine mka(a, min, max, diff) ! make array
    real(rk), allocatable :: a(:)
    real(rk) :: max, min, diff

    integer :: n, i

    n = (max - min)/diff

    allocate(a(n))

    a(1) = min
    do i = 2, n, 1
       a(i) = a(i - 1) + diff
    end do
  end subroutine mka

  subroutine column(r, c)
    integer :: c
    real(rk) :: r
    character :: ch
    character(3000) :: cch
    logical :: ws

    integer :: i, n, cc

    i = 1
    n = 3000
    cc = 0
    ws = .false.
    cch = ''

    do while (i < n)
       ch = line(i:i)
       !print *, trim(ch), trim(cch)

       if (ch == ' ' .or. ch == char(9) .or. ch == char(10)) then
          ws = .true.
       else
          if (ws .eqv. .true.) then
             ws = .false.
             cc = cc + 1

             if (cc == c) then
                exit;
             else
                cch = ''
             end if
          end if

          cch = trim(cch) // ch
       end if

       i = i + 1
    end do

    read(cch, *) r
  end subroutine column

  function parser(p, c) ! create array from file
    character(32) :: p
    integer :: c
    real(rk), allocatable :: parser(:)

    integer :: i, n, io
    real(rk), allocatable :: a(:)
    real(rk) :: s

    open(8, file = p, action = 'read')

    allocate(a(0))

    n = 2
    i = 1
    do while (i < n)
       read(8, '(A)', iostat = io) line
       call column(s, c)

       if (io == 0) then
          n = n + 1
          call append(a, s)
       end if
       i = i + 1
    end do

    allocate(parser(n))

    close(8)

    print *, a

    parser = a

    deallocate(a)
  end function parser

  subroutine append(a, v) ! append real to array
    real(rk), allocatable :: a(:)
    real(rk) :: v

    real(rk), allocatable :: b(:)
    integer :: n

    n = size(a)

    allocate(b(n))

    b = a
    deallocate(a)
    allocate(a(n + 1))

    a(1:n) = b
    a(n + 1) = v
  end subroutine append
end module plotutil

program main
  use plotutil
  use ogpf ! gnuplot

  implicit none

  character(32) :: s
  integer :: n

  real(rk), allocatable :: lum(:)
  real(rk), allocatable :: efftemp(:)
  real(rk), allocatable :: temp(:)
  real(rk), allocatable :: a(:)
  real(rk), allocatable :: b(:)
  real(rk), allocatable :: time(:)

  real(rk) :: s1 = 1.0140
  real(rk) :: a1 = 8.1774e-5
  real(rk) :: b1 = 1.7063e-9
  real(rk) :: c1 = -4.3241e-12
  real(rk) :: d1 = -6.6462e-16

  real(rk) :: s2 = 0.3438
  real(rk) :: a2 = 5.8942e-5
  real(rk) :: b2 = 1.6558e-9
  real(rk) :: c2 = -3.0045e-12
  real(rk) :: d2 = -5.2983e-16

  real(rk) :: solum = 3.828e26 ! W

  type(gpf):: gp

  call get_command_argument(1, s)

  print *, "Reading data..."

  time = parser(s, 4)

  efftemp = parser(s, 39) ! read data from file in command line arg

  lum = parser(s, 40) ! read data from file in command line arg

!  n = size()

  efftemp = 10**(efftemp) ! log -> lin
  lum = 10**(lum)

  temp = efftemp - 5780.0 ! t

  a = s1 + a1 * temp + b1 * temp**2 + c1 * temp**3 + d1 * temp**4 ! seff
  b = s2 + a2 * temp + b2 * temp**2 + c2 * temp**3 + d2 * temp**4

  a = (lum / a)**0.5 ! d
  b = (lum / b)**0.5

  !call mka(t, 0.0_8, real(n, kind = 8), 1.0_8)

  ! Annotation: set title, xlabel, ylabel
  call gp%title('Elinkelpoinen vyöhyke')
  call gp%xlabel('t')
  call gp%ylabel('d')
  call gp%options('set key top left; set style data linespoints')
  !Call Plot to draw a vector against a vector of data
  call gp%plot(x1=time,y1=a,ls1='title "moist"',x2=time,y2=b,ls2='title "max"',&
       x3=time,y3=a*0.95,ls3='title "moist -5%"',x4=time,y4=b*1.05,ls4='title "max +5%"')

!  call plot(t, v)
end program main