character infile*65,outfile*65,code(34)*1,line(105)*1,ans*1
      character deg*3,dg(3)*1,hyph(101)*1
      logical exists
      dimension num(98,98),jdeg(98)
      equivalence (deg,dg(1))
      data n,lsd,lnd,led,lwd,lesss,l,lesse,m,moren,morew,npt/12*0/
      data fs,fn,fe,fw/4*0.0/,hyph/101*'-'/
      data code/' ','1','2','3','4','5','6','7','8','9','A','B','C',
     1 'D','E','F','G','H','J','K','L','M','N','P','Q','R','S','T','U',
     2 'V','W','X','Y','Z'/
      print 100
  100 format(' INVENTPF.  Plouff, 4-94. Program to show regional grav',
     1 'ity data coverage by',/,'  printing the number of data points',
     2 ' in 2.5-minute cells.')
      print 101
  101 format (' TYPE the name of your data file (plouff format):')
      read (5,565) infile
  565 format (a65)
      inquire (file=infile,exist=exists)
      if (.not. exists) then
        print 102
  102   format (' Your file (<66 digits) was not found. Try again.')
        print 101
        read (5,565) infile
        inquire (file=infile,exist=exists)
        if (.not. exists) then
          print 103
  103     format (' **STOP.  Your file again was not found.')
          stop
        end if
      end if
      print 104
  104 format (' TYPE a name for a print file:')
      read (5,565) outfile
      inquire (file=outfile,exist=exists)
      if (exists) then
        print 105
  105   format (' Your print file already exists.  Do you want to ',
     1   'overwrite it (y,n)?')
        read (5,501) ans
  501   format (a1)
        if (ans .eq. 'y' .or. ans .eq. 'Y') go to 1
        print 104
        read (5,565) outfile
        inquire (file=outfile,exist=exists)
        if (exists) then
          print 106
  106     format (' **STOP.  That file also exists.')
          stop
        end if
      end if
    1 open (10,file=infile,status='old',form='formatted',blank='zero')
      open (15,file=outfile,status='unknown',form='formatted')
      write (15,100)
      write (15,600) infile
  600 format (' Input file: ',a65)
c Read data to find geographic range
    2 n=n+1
c Read geographic coordinates in plouff format, to 0.01 minute
      read (10,800,end=9,err=77) ltd,ltm,lgd,lgm
  800 format (bz,8x,i3,3i4)
      lt=ltm+6000*ltd
c Longitudes can be negative
      ln=lgm+6000*iabs(lgd)
      if (n .eq. 1) then
        minlt=lt
        minln=ln
        maxlt=lt
        maxln=ln
      else
        if (lt .lt. minlt) then
          minlt=lt
        else
          if (lt .gt. maxlt) maxlt=lt
        end if
        if (ln .lt. minln) then
          minln=ln
        else
          if (ln .gt. maxln) maxln=ln
        end if
      end if
      go to 2
   77 print 177, n
      write (15,177) n
  177 format (' **STOP. Format error in data file line',i5)
      go to 999
    9 n=n-1
      call conv (minlt,lsd,fs,lesss,m)
      call conv (maxlt,lnd,fn,l,moren)
      call conv (minln,led,fe,lesse,m)
      call conv (maxln,lwd,fw,l,morew)
      ncol=morew-lesse
      nrow=moren-lesss
c Print data limits.  Verify that array sizes are not exceeded.
      print 107, n,lsd,fs,lnd,fn,led,fe,lwd,fw
      write (15,107) n,lsd,fs,lnd,fn,led,fe,lwd,fw
  107 format (i7,' stations are in the input file.',/,'   Latitudes ',
     1 'range from',i4,'D',f6.2,'M to',i4,'D',f6.2,'M',/,'  Longitud',
     2 'es range from',i4,'D',f6.2,'M to',i4,'D',f6.2,'M')
c Geographic coordinates for 2.5-minute edges
      call conv25 (lesss,lsd,fs,moren,lnd,fn)
      call conv25 (lesse,led,fe,morew,lwd,fw)
c Test that array limits are not exceeded
      if (ncol .gt. 98) then
    3   print 108
  108   format (' Cell longitude range greater than the 4 degrees ',
     1   'is not permitted.',/,'  Type new map edges expressed to ',
     2   'the nearest 2.5 minutes.')
    4   print 109
  109   format (' TYPE the east edge coordinates (integer degrees, ',
     1   'decimal minutes):')
        read (5,*,err=4) led,fe
    5   print 110
  110   format (' TYPE the west edge coordinates (integer degrees, ',
     1   'decimal minutes):')
        read (5,*,err=5) lwd,fw
        write (15,601) led,fe,lwd,fw
  601   format (' Reduced cell borders are',i4,'D',f6.2,'M to',i4,
     1   'D',f6.2,'M longitude.')
        minln=iabs(led)+(100.0*fe+0.5)        
        maxln=iabs(lwd)+(100.0*fw+0.5)
        call conv (minln,led,fe,lesse,m)
        call conv (maxln,lwd,fw,l,morew)
        ncol=morew-lesse
        if (ncol .gt. 98) go to 3
      end if
      if (nrow .gt. 98) then
    6   print 111
  111   format (' Cell latitude range greater than the 4 degrees ',
     1   'is not permitted.',/,'  Type new map edges expressed to ',
     2   'the nearest 2.5 minutes.')
    7   print 112
  112   format (' TYPE the south edge coordinates (integer degrees,',
     1   ' decimal minutes):')
        read (5,*,err=7) lsd,fs
    8   print 113
  113   format (' TYPE the north edge coordinates (integer degrees,',
     1   ' decimal minutes):')
        read (5,*,err=8) lnd,fn
        write (15,602) lsd,fs,lnd,fn
  602   format (' Reduced cell borders are',i4,'D',f6.2,'M to',i4,
     1   'D',f6.2,'M longitude.')
        minlt=iabs(lsd)+(100.0*fs+0.5)        
        maxlt=iabs(lnd)+(100.0*fn+0.5)
        call conv (minlt,lsd,fs,lesss,m)
        call conv (maxlt,lnd,fn,l,moren)
        nrow=moren-lesss
        if (nrow .gt. 98) go to 6
      end if
      write (15,603) lsd,fs,lnd,fn,led,fe,lwd,fw
  603 format (' Map covers',i3,'D',f5.1,'M to',i3,'D',f5.1,'M lat and',
     1 i4,'D',f5.1,'M to',i4,'D',f5.1,'M long',/,' Numbers and let',
     2 'ters are the number of data points in a cell (no I, O; Z>32).')
c Initialize 2D array
      do 11 i=1,nrow
      do 10 j=1,ncol
   10 num(i,j)=0
   11 continue
c Add 1 to cell location, to convert to array location inside loop
      lesssm=lesss-1
      lessem=lesse-1
c Express map limits in 0.01-minute units
      lime=250*lesse
      lims=250*lesss
      limn=250*moren
      limw=250*morew
      rewind 10
      do 12 npts=1,n
      read (10,800) ltd,ltm,lgd,lgm 
      lt=ltm+6000*ltd
      ln=lgm+6000*iabs(lgd)
c Include south and east edges of map
      if (ln .lt. lime .or. ln .ge. limw .or. lt .lt. lims .or. lt .ge.
     1 limn) go to 12
      npt=npt+1
c Array relative to the SE corner. Truncation includes S, E cell edges.
      lt=(lt/250)-lesssm
      ln=(ln/250)-lessem
      num(lt,ln)=num(lt,ln)+1
   12 continue
c Position longitude integer degree locations in array
      ndeg=0
      do 13 j=1,ncol
c Number of 2.5-minute cells from Greenwich
      jcol=lessem+j
c Convert to degrees (truncated)
      jdegree=jcol/24
      if ((24*jdegree) .eq. jcol) then
        ndeg=ndeg+1
        jdeg(j)=jdegree
      else
        jdeg(j)=999
      end if
   13 continue
c Does west edge have integer degrees? yes=1=left
      jdegw=morew/24
      if ((24*jdegw) .eq. morew) then
        left=1
        ndeg=ndeg+1
      else
        left=0
      end if
c Total number of characters in each row
      limit=ncol+ndeg
c Mark longitude integer degrees at top
      limitop=limit
c Allow for 1 to 2 degree digit to extend beyond east edge
      if (jdeg(ncol) .ne. 999) limitop=limitop+2
      if (jdeg(ncol-1) .ne. 999) limitop=limitop+1
      do 14 j=1,limitop
   14 line(j)=' '
      jcol=0
      if (left .ne. 0) then
c Convert longitude degrees to alphameric characters
        call conval (jdegw,dg)
        line(1)=dg(1)
        line(2)=dg(2)
        line(3)=dg(3)
      end if
      do 15 j=1,ncol
c Reverse fill line from west to east (left to right)
      jr=ncol-j+1
      jcol=jcol+1
      if (jdeg(jr) .ne. 999) then
c Insert longitude degrees
        call conval (jdeg(jr),dg)
        line(jcol)=dg(1)
        jcol=jcol+1
        line(jcol)=dg(2)
        line(jcol+1)=dg(3)
      end if
   15 continue
      write (15,604) (line(j),j=1,limitop)
  604 format (1x,105a1)
      limitop=limit-11
c Mark north edge latitude if integer degrees
      jdegree=moren/24
      if ((24*jdegree) .eq. moren) write (15,605) jdegree,
     1  (hyph(j),j=1,limitop)
  605   format (i3,' degrees ',92a1)
c Loop latitudes, starting at north edge (reverse array)
      do 17 i=1,nrow
      ir=nrow-i+1
      jcol=0
      do 16 j=1,ncol
c Reverse array to start at west edge
      jr=ncol-j+1
      jcol=jcol+1
c Number of data points in cell plus one for array
      m=num(ir,jr)+1
      if (m .gt. 33) m=34
      line(jcol)=code(m)
      if (jdeg(jr) .eq. 999) go to 16
c Mark whole degree line east of its cell
      jcol=jcol+1
      line(jcol)='|'
   16 continue
      write (15,606) (line(k),k=1,limit)
  606 format (1x,103a1)
c Separate line for whole latitude degree line to south of its cell
      lt=ir+lesssm
      jdegree=lt/24
      if ((24*jdegree) .eq. lt) write (15,605) jdegree,
     1 (hyph(j),j=1,limitop)
   17 continue
      write (15,699) npt
      print 699, npt
  699 format (/,i7,' data points in the map.  Cells include ',
     1 'south and east edges.')
  999 close (10)
      close (15)
      stop
      end
      subroutine conv (lt,ltd,fm,less,more)
c Convert 0.01-minute geographic units (lt) to degrees/minutes and
c  2.5-minute cell edges (less/more in cell units, deg/min).
      ltd=lt/6000
      fm=0.01*(lt-6000*ltd)
      less=lt/250
      left=lt-250*less
      more=less
      if (left .ne. 0) more=less+1
      return
      end
      subroutine conv25 (l,ld,fl,m,md,fm)
      k=250*l
      ld=k/6000
      fl=0.01*(k-6000*ld)
      k=250*m
      md=k/6000
      fm=0.01*(k-6000*md)
      return
      end
      subroutine conval (num,dg)
c Converts absolute value of 3-digit numbers to 3 symbols
      character*1 dg(3),numb(10)
      data numb/'0','1','2','3','4','5','6','7','8','9'/
      n=iabs(num)
      n10=n/10
      n3=n-10*n10
      n1=n10/10
      n2=n10-10*n1
      dg(1)=numb(n1+1)
      dg(2)=numb(n2+1)
      dg(3)=numb(n3+1)
      return
      end