background map

[Dave Allured]

Enrico,

Grads map files must be written as three-byte records.  recl=3 is
correct, but your Fortran compiler must interpret the recl parameter as
single bytes, not "words".  This is a known standardization problem
between different compilers.  Your compiler may have an option to change
the meaning of the recl parameter.

Also, the use of the equivalence statement is platform dependent.  The
program requires the Big Endian number format.  Some compilers may have
a Little Endian option.  If that is not available, you can explicitly
reverse the byte ordering when writing.  For example, change this:

    write(20,rec=irec)(ilon(k),k=2,4)

to this:

    write(20,rec=irec) ilon(3), ilon(2), ilon(1)

and do the same for ilat.

I suspect that the record length is your main problem, but I don't have
enough information to be sure.  If you need a better diagnosis, then
please send a byte dump of the first part of your output file, using
this shell command:

    od -Ad -t u1 rupres | head

--Dave A.
CIRES Climate Diagnostics Center (CDC)
NOAA/ESRL, Physical Sciences Division (PSD)

enrico pisoni wrote:
> Dears
> I found a Fortran program on the GRADSUSR, kindly shared by Rupa Kumar
> Kolli (Indian Institute of Tropical Meteorology), that transform a file
> lon,lat of the following form:
>
>       2
>   6.000  41.000
>  17.000  47.000
>
> in a background map (i.e. named rupres), that it's possible to use in
> grads with:
> 'set mpdset rupres'
>
> But in grads I obtain the following error, invoking the 'set mpdset
> rupres':
>
> 'Map file format error: Invalid rec type 0 rec num 2'
>
> Someone can tell me which is the problem?
>
> I attach here Fortran Code I'm using
>
> -----------------------------------------------------------
>      PROGRAM GRADS_MAP
> c
> c     Program to prepare customized background map for GrADS from a
> c     text file containing line segments in (longitude,latitude)
> c     coordinates.
> c
> c     By Rupa Kumar Kolli, Indian Institute of Tropical Meteorology
> c        kolli at tropmet.ernet.in
> c        December 10, 1997
> c
> c     Input File Format:
> c     The input file can contain any number of line segments.
> c     Each segment starts with the number of pairs of coordinates (I8)
> c     The coordinates should be listed from the subsequent record,
> c         as longitude,latitude pair (2F8.3), one pair per line.
> c
> c     Example:
> c
> c           5
> c     209.830 -76.580
> c     211.500 -76.750
> c     212.830 -76.420
> c     211.830 -76.170
> c     209.830 -76.580
> c
> c     VERY IMPORTANT: NOTE THE FOLLOWING !!!
> c
> c     1. Longitudes should be in the range 0 to 360 (No negative
> c        values !).
> c
> c     2. Latitudes should be in the range -90 (South Pole) to 90
> c        (North Pole).
> c
> c     3. DON'T CROSS GREENWICH MERIDIAN (0.0/360.0) IN ANY SEGMENT.
> c        For lines across this meridian, split them into two parts,
> c        with a common point between them.  Use the value 360.0 for
> c        Greenwich meridian if you are moving towards it from west to
> c        east, and the value 0.0 if you are moving from east to west.
> c        That is, you cannot jump from 359.0 to 0.0; this will wrap
> c        the line across the whole globe !
> c
> c     4. The maximum number of points (pairs of coordinates) in a
> c        single segment is 255.
> c
> c     5. Don't use GrADS standard file names (lowres, mres and hires)
> c        for the output file (e.g., rupres). Preserve them for
> c        standard use.
> c
> c     6. After generating the GrADS map background file, copy it to a
> c        path where GrADS looks for its support files (generally the
> c        path defined in $GADDIR)
> c
> c     7. Finally, start GrADS and give the command
> c        'set mpdset filename':
> c        e.g.,
> c        set mpdset rupres
> c
> c     8. This program was developed on a HP 9000/720 machine running
> c        HP-UX 9.0 operating system.  It could work on other
> c        platforms as well, but I have not tested.
> c
> c     ENJOY HAVING YOUR OWN LINES !!!
>
>      byte ihead(3),ilon(4),ilat(4)
>      integer lon,lat
>      equivalence (lon,ilon), (lat,ilat)
> c
> c     Open the input and output files
> c
>      open(10,file='rupres.txt',status='old')
>      open(20,file='rupres',form='unformatted',access='direct',recl=3)
> c
> c     Initialize the first two header-record bytes: always 1
> c
>      ihead(1)=1
>      ihead(2)=1
>      irec=0
>      nsegs=0
> c
> c     Read the number of pairs in this segment.
> c     Leave the loop if end of file is reached.
> c
>   99 read(10,1,end=100)npairs
>    1 format(i8)
>      nsegs=nsegs+1
>      if (npairs.gt.255) then
>       print *,'Number of pairs in Segment No.',nsegs,' exceeds 255 !'
>       stop 'Job Aborted'
>      endif
> c
> c     Write the header record (3 bytes)
> c     The third byte indicates the number of pairs
> c     of values following (maximum=255)
> c
>      ihead(3)=npairs
>      irec=irec+1
>      write(20,rec=irec)ihead
> c
> c     Loop over the current segment
> c
>      do 10 i=1,npairs
> c
> c     Read one pair of coordinates (Longitude,Latitude)
> c
>        read(10,2)alon,alat
>    2   format(2f8.3)
> c
> c     Convert the values to integers, in GrADS readable format
> c
>        lon=int(alon*10000.)+0.5
>        lat=int((alat+90.)*10000.)+0.5
> c
> c     Write the last three-bytes of lon and lat on the output file
> c
>        irec=irec+1
>        write(20,rec=irec)(ilon(k),k=2,4)
>        irec=irec+1
>        write(20,rec=irec)(ilat(k),k=2,4)
>   10 continue
> c
> c     Go back to read next line segment
> c
>      goto 99
>  100 print *,'Total Number of Segments Processed=',nsegs
>      close (10)
>      close (20)
>      stop 'Normal Termination'
>      end
> ----------------------------------------------------------