background map

[Davide Sacchetti]

supposing input format to be correct, the problem could be on the output
file: in the fortran open statement, the RECL keyword behaves
differently upon different fortran compilers ...
here you have RECL=3
try with RECL=12

bye bye
Davide


On Fri, 2006-07-28 at 10:52 +0200, 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$FGPR
> 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$FGPR
> 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
> ----------------------------------------------------------
--
Sacchetti Davide
ARPAL UO3 Centro Meteo Idrologico Regione Liguria - Dir. Scientifica
V.le Brigare Partigiane 2 16121 Genova (I)
tel: +39 010 6437535
mail: davide.sacchetti at arpal.org   web: www.meteoliguria.it