Gravsoft code - Fortran 95

Hello,

I have a simple program that reads a datafile in scanline format and writes out an XYZ file (longitude, latitude and value). Everything is set fine and the input and output files are read fine. I have modified existing code to read the Gravsoft data format.

However, I would like to make the dimension statement generic, so that the user can input a new value if the input file has a larger dimension (eg more columns).

If I put Implicit none in, then I have to define nn, ne, i, j, and row as integer. The code compiles but a run-time error is generated, so this is left out to make it work.

Any pointers how to streamline the code and make it more efficient would be helpful.

What I am looking to achieve is a generic Gravsoft reader, e.g.:

Enter input data: Enter output data: Enter row dimension:

Program Gravsoft
dimension row(4320)
real :: rlat,rlat1,rlat2,rlon,rlon1,rlon2,dlat,dlon
character(len=128) :: infile, outfile
print *, 'Enter input data file:'
read (*,'(A40)') infile
open(10, file=infile, status='old')
print *, 'Enter output data file:'
read (*,'(A40)') outfile
open(11, file=outfile, status='new')
read(10,*) rlat1,rlat2,rlon1,rlon2,dlat,dlon
nn = (rlat2-rlat1)/dlat+1.5
ne = (rlon2-rlon1)/dlon+1.5
rlat = rlat2
do i = 1, nn
  read(10,*) (row(j),j=1,ne)
  rlon = rlon1
  do j = 1,ne
    if (rlon .gt. 360) rlon = rlon1
    write(11,*) rlat, rlon, row(j)
    rlon = rlon + dlon
  enddo
  rlat = rlat - dlat
enddo
close (10)
close (11)
end program

Thanks

Lester

For a start, you could make the dimension a lot more. You are using REAL, so your ROW array takes 43204 bytes = 17280 (or 34560 if you used REAL8) - and you aren't going to find a modern computer without at least 1000,000,000 free for a program like this ...

More to the point, why is ROW an array anyway? A single variable would do the job and be faster (although marginally so with an array of this size).

If you are using 4-byte REALs, then incrementing/decrementing RLAT and RLON will eventually give you roundoff, and you might be better doing calculations with i and j.

Do you mean something like:

Program Gravsoft
   
   implicit none
   
   ! Local Variables
   integer i, j, nn, ne
   real rlat,rlat1,rlat2,rlon,rlon1,rlon2,dlat,dlon
   real, allocatable :: row(:)
   character(len=128) :: infile, outfile
   
   ! Read file names from keyboard
   print *, 'Enter input data file:'
   read (*,'(A)') infile
   print *, 'Enter output data file:'
   read (*,'(A)') outfile
   
   ! Open input and output files
   open(10, file=infile, status='old', action='read', position='rewind')
   open(11, file=outfile, status='new', action='write')
   
   ! Read grid data and calculate number of rows (nn) and columns (ne)
   read(10,*) rlat1,rlat2,rlon1,rlon2,dlat,dlon
   nn = int((rlat2-rlat1)/dlat+1.5)
   ne = int((rlon2-rlon1)/dlon+1.5)
   
   ! Create space for data on each row
   allocate(row(ne))
   
   rlat = rlat2
   do i = 1, nn
      read(10,*) (row(j),j=1,ne)
      rlon = rlon1
      do j = 1, ne
         if (rlon .gt. 360.0) rlon = rlon1
         write(11,*) rlat, rlon, row(j)
         rlon = rlon + dlon
      enddo
      rlat = rlat - dlat
   enddo
   
   ! Destroy space for data (not essential in F95)
   deallocate(row)
   
   ! Close files
   close (10)
   close (11)
end program

Agree with above though; you should really calculate rlon and rlat from the i,j indices. Aside for correcting for round-off, this also opens up the possibility of parallelising the loops for a compiler with such a capability.

Hi David,

Yes the code you show is what I was looking for, thanks for the corrections and the pointers - always helpful. The orginal code came from http://earth-info.nga.mil/GandG/wgs84/agp/readme_new.html and definitely not very elegant.

The code works fine - and no need to explicitly cite the row dimension. The only other change I made was to set rlat,rlon etc to real*8 (when scanning the header line in the data file) as this is positional data at fine resolution.

Cheers Lester

(http://earth-info.nga.mil/GandG/wgs84/agp/readme_new.html) and definitely not very elegant.

The code works fine - and no need to explicitly cite the row dimension. The only other change I made was to set rlat,rlon etc to real*8 (when scanning the header line in the data file) as this is positional data at fine resolution.

Cheers Lester

[quote:fa16a87e3f="davidb"]Do you mean something like:

Program Gravsoft
   
   implicit none
   
   ! Local Variables
   integer i, j, nn, ne
   real rlat,rlat1,rlat2,rlon,rlon1,rlon2,dlat,dlon
   real, allocatable :: row(:)
   character(len=128) :: infile, outfile
   
   ! Read file names from keyboard
   print *, 'Enter input data file:'
   read (*,'(A)') infile
   print *, 'Enter output data file:'
   read (*,'(A)') outfile
   
   ! Open input and output files
   open(10, file=infile, status='old', action='read', position='rewind')
   open(11, file=outfile, status='new', action='write')
   
   ! Read grid data and calculate number of rows (nn) and columns (ne)
   read(10,*) rlat1,rlat2,rlon1,rlon2,dlat,dlon
   nn = int((rlat2-rlat1)/dlat+1.5)
   ne = int((rlon2-rlon1)/dlon+1.5)
   
   ! Create space for data on each row
   allocate(row(ne))
   
   rlat = rlat2
   do i = 1, nn
      read(10,*) (row(j),j=1,ne)
      rlon = rlon1
      do j = 1, ne
         if (rlon .gt. 360.0) rlon = rlon1
         write(11,*) rlat, rlon, row(j)
         rlon = rlon + dlon
      enddo
      rlat = rlat - dlat
   enddo
   
   ! Destroy space for data (not essential in F95)
   deallocate(row)
   
   ! Close files
   close (10)
   close (11)
end program

Agree with above though; you should really calculate rlon and rlat from the i,j indices. Aside for correcting for round-off, this also opens up the possibility of parallelising the loops for a compiler with such a capability.