!
! To compile on Alpha machine:
!    f90 -free -cpp -convert big_endian read_veg.F
! To complie on Linux:
!    pgf90 -Mfreeform -DRECLENBYPTE -byteswapio read_veg.F
!
program read_vegetation

! Open direct access file VEG-USGS.30s, and read it
!      Note 30 sec vegetation data starts at 
!           latitude 90. and longitude -180.
!
!    - the file is in direct access format
!    - each direct access record contains data in one latitude circle
!      beginning at -180 degree longitude, and end at +180 degree
!    - the data is arranged to start at northernmost latitude (north pole),
!      and end at south pole

 implicit none

! declare variables

 character (len=25) :: vegfile
 character (len=1 ), allocatable, dimension(:) :: veg_char

 integer :: iunit = 10
 integer :: rec_len_lat = 360*120
 integer :: rec_len_lon = 180*120 
 integer :: rec_len, length, irec, jrec, nrec, ierr
 integer :: begin_lat_rec, end_lat_rec, begin_lon_rec, end_lon_rec
 real    :: begin_lat, end_lat, begin_lon, end_lon
 integer, allocatable, dimension(:) :: veg_int

! get the file name from input

 call getarg(1,vegfile)
 print *, 'reading file ', vegfile

! area of data to be printed
! longitudes are bounded between -180 W to +180 E

 begin_lat = 50.
 begin_lon = -130.
 end_lat   = 30.
 end_lon   = -70.

! record length for the data
!    each record has 360x120 data points

 rec_len = rec_len_lat
 length  = rec_len/4

#ifdef RECLENBYTE                                
! modify for machines like Cray, Sun, HP, IBM and Linux
 length = length*4                      
#endif                                         

! open 30 sec vegetation data file: 

 open (iunit,file=vegfile,access='DIRECT',recl=length,status='OLD')

! begin_record and end_record may be calculated as follows:
!
 begin_lat_rec = nint(90.-begin_lat)*120 + 1
   end_lat_rec = nint(90.-  end_lat)*120 + 1
 begin_lon_rec = nint(begin_lon-(-180.))*120 + 1
   end_lon_rec = nint(end_lon  -(-180.))*120 + 1

 print *, 'begin_lat_rec and end_lat_rec = ', begin_lat_rec, end_lat_rec
 print *, 'begin_lon_rec and end_lon_rec = ', begin_lon_rec, end_lon_rec

! nrec: from north to south (max 21600)
! irec: from west to east   (max 43200)

 allocate (veg_char(rec_len))
 allocate (veg_int (rec_len))

!do nrec = begin_lat_rec, end_lat_rec
 do nrec = 1, rec_len_lon

    read (iunit, rec=nrec, iostat=ierr) veg_char

    if (ierr /= 0) then
       print *, 'read error on record nrec = ', nrec
       stop 'exit read error'
    endif

    veg_int = ichar (veg_char)
    if ( nrec >= begin_lat_rec .and. nrec <= end_lat_rec) then
    if (mod(nrec,120) == 0 ) print *, (veg_int(irec), irec=begin_lon_rec,end_lon_rec,120)
    end if

 end do

 deallocate (veg_char)
 deallocate (veg_int )

 stop
end program read_vegetation