Why are incorrect values printed outside of the Fortran loop?

Question

I read 3 arrays (X_halo, Y_halo, Z_halo) from an input file. If I write the values loaded within the DO cycle I can print them correctly, but out of the loop no. And if I try to print the values of X_halo(i) incorrect values are printed.

How can I do it?

program columns
  IMPLICIT NONE

  INTEGER,SAVE :: lun
  INTEGER, PARAMETER :: ARRAYLEN=1044  
  CHARACTER(120) :: filename
  DOUBLE PRECISION, DIMENSION (ARRAYLEN) :: X_halo, Y_halo, Z_halo
  INTEGER :: i, istat

  lun=1
  filename = 'xyz.dat'

  OPEN (UNIT=10, FILE=filename, STATUS='old', ACTION='read', IOSTAT=istat)

    DO i=1,ARRAYLEN
     READ (10, *, iostat=istat) X_halo(i), Y_halo(i), Z_halo(i)
!    print*, 'id:', i, 'X= ', X_halo(i), 'Y= ', Y_halo(i), 'Z= ', Z_halo(i)
    END DO

  DO i=1,ARRAYLEN
  print*, X_halo(i)
  END DO

  CLOSE (10)

end program columns

For example the first lines of the first column of the file xyz.dat are:

But the code returns:

6.9532597489392050E-310
2.2310395176993305E-314
6.9532121310250636E-310
6.9532238136146167E-310

Second problem: If you want to have ARRAYLEN as free variable to use the program for all the files in put, how can I do it?

Pleasemake your titles like real short questions. And do not add the tag. See [ask] and https://meta.stackexchange.com/questions/19190/should-questions-include-tags-in-their-titles — Vladimir F Героям слава, Aug 10 '17 at 18:16
If you have two problems, consider aking two separate questions. (I have not read your current problems yet.) — Vladimir F Героям слава, Aug 10 '17 at 18:17
Also, you add the [tag:fortran90] tag to your question. Are you sure you want to limit yourself to this very old version? There is Fortran 95, 2003, 2008 and 2015. Fortran 90 is really old, do you require sticking to it? If not, just use [tag:fortran] only. — Vladimir F Героям слава, Aug 10 '17 at 18:19
Which numbers do you read? What is printed? It is always good to describe exactly what is wrong. Regarding your other point, read about allocatable arrays. — Vladimir F Героям слава, Aug 10 '17 at 18:22
There are several problems with your program. You do not check `istat`, you close unit `` which has not been opened, but please explain what is wrong with the output. — Vladimir F Героям слава, Aug 10 '17 at 18:28
You still have the problem of not checking the value of `istat`. You **have to** check that! Or delete the `iostat=istat`. If the file os too short, you have a big problem. Try setting `ARRAYLEN=4`. — Vladimir F Героям слава, Aug 10 '17 at 19:49
See these related questions and answers https://stackoverflow.com/questions/32682309/number-of-lines-of-a-text-file https://stackoverflow.com/questions/44371378/count-lines-within-a-fortran-loop https://stackoverflow.com/questions/30692424/how-to-read-number-of-lines-in-fortran-90-from-a-text-file — Vladimir F Героям слава, Aug 10 '17 at 19:56
I tried to delete stat and now it seems to work, but why does it create problems? — Alessandro Peca, Aug 10 '17 at 20:36
Well in that case you should have received an error message. Did you? We have to see the **complete** output. It looks you are not showing everything. — Vladimir F Героям слава, Aug 10 '17 at 20:41
FYI, the variable `lun` is unused in this example. Note that typically this would stand for "logical unit". Since it is defined you may as well use it as in `open(unit=lun, ...)`. — Matt P, Aug 11 '17 at 15:23

Steve · Accepted Answer · 2017-08-10T23:35:47.073

Not sure why you're receiving an error due to istat or why your values aren't printing correctly - they're fine for me, will need a minimal complete verifiable example with error messages to get feedback on that.

With regards to reading fewer lines than ARRAYLEN, as discussed in the comments in your previous question and outlined in the comments above, you want to check if istat < 0 then break from the loop, with

IF ( istat < 0 ) THEN
  EXIT
END IF

One way of having X_halo to be the length of lines read (I doubt this is the best way) is to save the entries to temporary vectors, check how many entries were read in with i, ALLOCATE your halo vectors, then copy from the temporary vectors.

PROGRAM columns
  IMPLICIT NONE

  INTEGER,SAVE :: lun
  INTEGER, PARAMETER :: ARRAYLEN=1044
  CHARACTER(120) :: filename
  ! Read data into tempory arrays of length up to ARRAYLEN
  DOUBLE PRECISION, DIMENSION (ARRAYLEN) :: X_tmp, Y_tmp, Z_tmp
  ! And these will be where we put the data when we know how many lines
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION (:) :: X_halo, Y_halo, Z_halo
  INTEGER :: i, istat

  lun = 1
  filename = 'xyz.dat'

  OPEN (UNIT=10, FILE=filename, STATUS='old', ACTION='read', IOSTAT=istat)

  DO i=1,ARRAYLEN
    ! Read data in to the temporary vectors
    READ (10, *, IOSTAT=istat) X_tmp(i), Y_tmp(i), Z_tmp(i)
    IF ( istat < 0 ) THEN
      ! Then last READ failed with an EOF
      EXIT
    END IF
  END DO

  ! Last READ returned EOF, so i-1 lines were read
  i = i-1

  ALLOCATE(X_halo(i))
  ALLOCATE(Y_halo(i))
  ALLOCATE(Z_halo(i))
  ! Copy i entries in to the allocated vectors
  X_halo = X_tmp(1:i)
  Y_halo = Y_tmp(1:i)
  Z_halo = Z_tmp(1:i)

  PRINT*, X_halo

  CLOSE (10)

END PROGRAM columns

Running this with file xyz.dat as

281.0788 1.3 111
189.8768 2.3 222
669.2193 3.3 333
720.7653 4.3 444

returns

   281.07880000000000        189.87680000000000        669.21929999999998        720.76530000000002

Another approach would be to READ through the file not saving the output just to get the number of lines i, then ALLOCATE your vectors, REWIND and READ in each line. See How to read number of lines in Fortran 90 from a text file?

Thank you @Steve. I made a little change in `y = i-1` and `ALLOCATE(X_halo(y))` `X_halo = X_tmp(1:y)` so i can do `DO i=1,y PRINT*, X_halo(i) ENDDO` to print in column. — Alessandro Peca, Aug 11 '17 at 06:42
If you go for the "read to count lines, rewind, read to save data" approach you don't need `ARRAYLEN` for max lines (`DO` loop without condition) or the `temp` vectors at all. — Steve, Aug 11 '17 at 07:13

Why are incorrect values printed outside of the Fortran loop?

1 Answers1