|
forums.silverfrost.com Welcome to the Silverfrost forums
|
View previous topic :: View next topic |
Author |
Message |
arctica
Joined: 10 Sep 2006 Posts: 105 Location: United Kingdom
|
Posted: Tue Oct 21, 2014 9:08 am Post subject: 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:
Code: | 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 |
|
Back to top |
|
|
LitusSaxonicum
Joined: 23 Aug 2005 Posts: 2388 Location: Yateley, Hants, UK
|
Posted: Tue Oct 21, 2014 5:11 pm Post subject: |
|
|
For a start, you could make the dimension a lot more. You are using REAL, so your ROW array takes 4320*4 bytes = 17280 (or 34560 if you used REAL*8) - 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. |
|
Back to top |
|
|
davidb
Joined: 17 Jul 2009 Posts: 560 Location: UK
|
Posted: Tue Oct 21, 2014 6:17 pm Post subject: |
|
|
Do you mean something like:
Code: |
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. _________________ Programmer in: Fortran 77/95/2003/2008, C, C++ (& OpenMP), java, Python, Perl |
|
Back to top |
|
|
arctica
Joined: 10 Sep 2006 Posts: 105 Location: United Kingdom
|
Posted: Wed Oct 22, 2014 9:44 am Post subject: Re: |
|
|
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
davidb wrote: | Do you mean something like:
Code: |
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. |
|
|
Back to top |
|
|
|
|
You cannot post new topics in this forum You cannot reply to topics in this forum You cannot edit your posts in this forum You cannot delete your posts in this forum You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|