Issues with Sequential Call of Fortran Subroutines

Question

I have a Fortran script that calls multiple subroutines from an internal module (see below).

The subroutines are called sequentially.

I noticed that the first two calls are executed normally, but the last subroutine is not executed.

I have tried converting the subroutines into functions, but run into the same issue.

Will appreciate ideas on how to solve this issue.

module tmod

  contains
  subroutine readata(fname,tm,hd,td,nrec)
    ! read input data and returns arrays
    implicit none
    integer :: io,i
    integer,intent(in) :: nrec
    character(len=40) :: fname
    real(kind=8),allocatable,intent(out) :: tm(:),hd(:),td(:)
    
    allocate(tm(nrec-1))
    allocate(hd(nrec-1))
    allocate(td(nrec-1))
    
    open(10,file=fname,status='old',action='read')
      read(10,*)
      do i=1,nrec-1
        read(10,*) tm(i),hd(i),td(i)
      end do
    close(10)

    return
  end subroutine readata

  subroutine fd(nrec,t_arr,v_arr,derv_arr)
    implicit none
    integer,intent(in) :: nrec
    integer :: i
    real(kind=8),allocatable,intent(out) :: derv_arr(:)
    real(kind=8),dimension(:),intent(in) :: t_arr, v_arr
  
    allocate(derv_arr(nrec-2))   

    do i=2,nrec-1
      derv_arr(i)= (v_arr(i+1)-v_arr(i-1))/(t_arr(i+1)-t_arr(i-1))
    end do
    
    return
  end subroutine fd

  subroutine sd(nrec,t_arr,v_arr,derv2_arr)
    implicit none
    integer,intent(in) :: nrec
    integer :: i
    real(kind=8),dimension(:),intent(in) :: t_arr,v_arr
    real(kind=8),allocatable,intent(out) :: derv2_arr(:)

    allocate(derv2_arr(nrec-2))

    do i=2,nrec-1
      derv2_arr(i)=(v_arr(i+1)-2*v_arr(i)+v_arr(i-1))
    end do

    return

  end subroutine sd
end module tmod

program main
  use tmod
  implicit none
  character(len=40) :: filename
  real(kind=8),allocatable :: tg(:), hg(:), tig(:),pd2(:), sed(:), pd1(:)

  integer :: nrg

  ! input parameters
  filename="test_data.dat"
  nrg=5
  
  ! calling subroutines
  call readata(filename, tg, hg, tig, nrg)
  write(*,*) hg

  call fd(nrg,tg,tig,sed)
  write(*,*) "First"
  write(*,*) sed

  call fd(nrg,tg,hg,pd1)
  write(*,*) "First 1"
  write(*,*) pd1
 
  call sd(nrg,tg,hg,pd2)
  write(*,*) "Second"
  write(*,*) pd2
  
end program maintype

Below is the content of the input file "test_data.dat"

t   h   ti
0   12.113  1.1e-3
15  12.123  1.12e-3
30  12.156  1.11e-3
45  12.134  1.15e-3

Answer 1

If you'd have compiled your application with array bounds checking enabled, you'd have found your problem quickly:

$ gfortran -O2 -fcheck=all -g callbug.f90 
$ ./a.out 
   12.113000000000000        12.122999999999999        12.156000000000001        12.134000000000000     
At line 36 of file callbug.f90
Fortran runtime error: Index '5' of dimension 1 of array 'v_arr' above upper bound of 4

Your fundamental problem seems to be that your nrec variable in the subroutines is not the number of records you have, but the number of lines in your data files. The actual number of records being one less due to the header line in the file, but then since you're using the i+1 index in the loop body you must stop the iteration already at nrec-2. One way to avoid these kinds of mistakes, which are all too easy to make, apart from enabling the bounds checking, is to use the Fortran size intrinsic to get the size of the array you're iterating over instead of a separate nrec variable. So for instance, on lines 35 and 51 in your program, instead of

do i=2,nrec-1

write it as

do i=2,size(v_arr)-1

Similarly when allocating the size of the result array on lines 33 and 49, instead of

allocate(derv_arr(nrec-2))

use

allocate(derv_arr(size(v_arr)-2))

Now you can remove the unused nrec argument from the fd and sd subroutines.

Another issue is that as you're starting the iterations from i=2, you're leaving the first elements of the result arrays unused, so they're containing some nonsense values. Thus in your loops you should have derv_arr(i-1)=....

You can also avoid having to hardcode the number of lines in your data file by first counting the number of lines, rewriting readdata as follows:

  subroutine readata(fname,tm,hd,td)
    ! read input data and returns arrays
    use, intrinsic :: ISO_FORTRAN_ENV
    implicit none
    integer :: i, istat, nrec, u
    character(len=*) :: fname
    real(kind=real64),allocatable,intent(out) :: tm(:),hd(:),td(:)
      
    open(newunit=u,file=fname,status='old',action='read')
    read(u, *)
    nrec = 0
    do
      read(u, *, iostat=istat)
      if (istat /= 0) then
        exit
      end if
      nrec = nrec + 1
    end do
    allocate(tm(nrec))
    allocate(hd(nrec))
    allocate(td(nrec))

    rewind(u)
    read(u, *)
    do i=1,nrec
      read(u, *) tm(i), hd(i), td(i)
    end do
    close(u)
  end subroutine readata

Note a few other details:

• Use character(len=*) for character arguments to automatically get the size from the caller, avoiding yet another common source of running beyond the bounds of the variable.
• use, intrinsic :: ISO_FORTRAN_ENV and kind=real64 instead of kind=8 for real variables. You should do this everywhere, of course, not only in this subroutine.
• Use newunit= instead of hardcoding a unit number.