      PROGRAM GRADS_MAP
c
c     Program to prepare customized background map for GrADS from a text 
c     file containing line segments in (longitude,latitude) coordinates.
c
c     By Rupa Kumar Kolli, Indian Institute of Tropical Meteorology
c        kolli@tropmet.ernet.in
c        created:  December 10, 1997
c        modified: March 12, 1998
c        PC Fortran support added: January 26, 1999
c     Please do contact me for comments/feedback etc. on this program.
c
c     This program converts map files in ASCII into GrADS readable format.
c     The resolution of the map background, of course, depends on the 
c     basic data that you can get.  The format of the GrADS map data set
c     prepared by this program conforms to that used in 'lowres'.  While it
c     should work for any kind of resolution, it may be a bit slower than 
c     the standard 'hires' for huge files.
c
c     Input File Format:
c     The input file can contain any number of line segments.
c     Each segment starts with a header record containing the
c         the number of pairs of coordinates and associated map type (2I8)
c         (the map type is used to selectively display map lines using
c          the 'set mpt ...' command)
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       2
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 values !).
c
c     2. Latitudes should be in the range -90 (South Pole) to 90 (North Pole).
c
c     3. DON'T CROSS THE GREENWICH MERIDIAN (0.0 or 360.0) IN ANY SEGMENT.
c        For lines across this meridian, split them into two parts, with
c        a common point between them.  Use the value 360.0 for Greenwich
c        meridian if you are moving towards it from west to east, and the
c        value 0.0 if you are moving from east to west.  That is, you
c        cannot go from 359.0 to 0.0; this will wrap the line across the
c        whole globe !
c
c     4. The maximum number of points (pairs of coordinates) in a single
c        segment is 255.
c
c     5. Don't use GrADS standard file names (lowres, mres, hires, nam) for
c        the output file (e.g., rupres). Preserve them for standard use.
c
c     6. After generating the GrADS map background file, copy it to a path
c        where GrADS looks for its support files (generally the path defined
c        in $GADDIR)
c
c     7. Finally, start GrADS and give the command 'set mpdset filename':
c        e.g., set mpdset rupres.  Note that you can keep multiple map
c        data sets open in GrADS 1.7; thus, you can create separate map
c        data files and selectively use subsets of the files.  For this
c        purpose, use the command like 'set mpdset map1 map2', whereupon
c        GrADS uses both the files for drawing the map.
c
c     8. This program was initially developed on a HP 9000/720 machine with
c        HP-UX 9.05 operating system.  I have successfully tested it in dos/
c        Windows as well as in RedHat Linux (g77). It should work on other
c        platforms as well.  For using this program on a PC (dos/windows/
c        linux), a slight change is necessary, as indicated in the comments
c        preceding the output statements (marked by --->) in the program.
c
c     ENJOY HAVING YOUR OWN CUSTOMIZED GrADS MAP !!!

      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='mexico_s.asc',status='old')
      open(20,file='mexico_s',form='unformatted',
     &              access='direct',recl=3)
c
c     Initialize the first header-record byte, ihead(1): fixed to be 1
c     A value of 2 can be used to skip unwanted sections of the file,
c     (mres and hires were created like this), but this facility
c     is not implemented in this program.  So don't change this value
c     unless you want to implement skipping (could be fairly complicated !)
c
      ihead(1)=1
      irec=0
      nsegs=0
c
c     Read the number of pairs and map type for this segment.
c     Leave the loop if end of file is reached.
c
   99 read(10,1,end=100)npairs,ihead(2)
    1 format(2i8)
      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
      write(*,*)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
c        write(*,*)npairs,alon,alat
cjM        read(10,*)alon,alat
    2   format(2f11.5)
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     For use with PCs running DOS/Windows/Linux:
c     Write the first three-bytes of lon and lat on the output file,
c     in reverse order.
c     Comment out the statements meant for Workstations.
c
c     For use with HP-UX, IRIX, SOLARIS:
c     Write the last three-bytes of lon and lat on the output file.
c     Comment out the statements meant for PC.
c     
c     Note that ilon and ilat are byte equivalents of lon and lat.
c
        irec=irec+1
c --->  For PC (dos/windows/linux):
        write(20,rec=irec)(ilon(k),k=3,1,-1)
c --->  For Workstations (RISC):
c        write(20,rec=irec)(ilon(k),k=2,4)
        irec=irec+1
c --->  For PC (dos/windows/linux):
        write(20,rec=irec)(ilat(k),k=3,1,-1)
c --->  For Workstations (RISC):
c        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
