C     program shefplot.exe
      common/limits/ylatn,xlonw,ylats,xlone,xof,yof
      dimension iasf(13),ux(15000),uy(15000),iuv(15000)
      dimension          hx(500), hy(500),  kuv(500)
      dimension          px(500), py(500),  luv(500)
      dimension          ax(500), ay(500),  muv(500)
      integer divsor, roundr, hi, pi, ai
      character *62 infile
      character *60 title, fecha

      do i=1,13
         iasf(i)=1
      enddo
      call opngks
      call gsclip(0)
      call gsasf(iasf)

c---  read in the file name for the input file, variable infile
      read(5,25) infile
   25 format(a62)
      open(unit=11,file=infile,form='formatted',status='old')

C Use EZMAP and EZMAPA to create a background.

      call mappos (0.,1.,0.,1.)
      call maproj ('CE - stereo',60.,-99.,0.)
C     call maproj ('CE - stereo',60.,-100.,0.)
      ylatn= 50.
      xlonw= -131.
      ylats= 20.
      xlone= -74.
C     xlone= -76.
      xof  = .075
      yof  = .12
      call mapset ('CO - CORNERS',ylatn,xlonw,ylats,xlone)
      call mapsti ('EL - ELLIPTICAL BOUNDARY',0)
      call mapsti('C7',0)
      call mapsti('C6',1)
      call mapstc ('OU - OUTLINE DATASET','CO')
      call mapint
      call getset (xvpl,xvpr,yvpb,yvpt,xwdl,xwdr,ywdb,ywdt,lnlg)
C
C Outline the continents in black, put black contour lines over the
C color map, and put gray lines of latitude and longitude over the ocean
 
      call gsplci (1)
      ii=0
      hi=0
      pi=0
      ai=0
      read(11,40) title, fecha
   40 format(a60/a60)
      write(*,*) title
      write(*,*) fecha
      read(11,*)mark,divsor 
 	
c     determine the rounding factor by multiplying .5
c     times the divsor. if no divsor on input card, set to 1
 
      if (divsor.ne.1.and.divsor.ne.10.and.divsor.ne.100) divsor=1
      roundr = .5 * divsor
  101 continue
c
c     read and load array with the sites in the 48 mainland US states
c 
      read(11,*,end=500) y,x,iaz 
      if (y .lt. 20.0 .or. y .gt.  50.0) go to 201
      if (x .lt. 65.0 .or. x .gt. 127.0) go to 201
      iaz = (iaz + roundr) / divsor
      x=x*(-1.)
      if(iaz.lt.1999)then
         ii=ii+1
         ux(ii)=x
         uy(ii)=y
         iuv(ii)=iaz
c     write (*,444) ii,ux(ii),uy(ii),iuv(ii)
 444  format (i5,2f8.2,2i5)
      endif
      go to 101
c
c     read and load array with the sites in Hawaii  
c 
  201 continue
      if (y .lt. 18.5 .or. y .gt.  22.5) go to 301
      if (x .lt.154.0 .or. x .gt. 161.0) go to 301
      iaz = (iaz + roundr) / divsor
      x=x*(-1.)
      if(iaz.lt.1999)then
         hi=hi+1
         hx(hi)=x
         hy(hi)=y
         kuv(hi)=iaz
      endif
      go to 101
c
c     read and load array with the sites in Puerto Rico
c 
  301 continue
      if (y .lt. 17.5 .or. y .gt. 19.0) go to 401
      if (x .lt. 65.0 .or. x .gt. 67.5) go to 401
      iaz = (iaz + roundr) / divsor
      x=x*(-1.)
      if(iaz.lt.1999)then
         pi=pi+1
         px(pi)=x
         py(pi)=y
         luv(pi)=iaz
      endif
      go to 101
c
c     read and load array with the sites in Alaska     
c 
  401 continue
      if (y .lt. 54.0 .or. y .gt. 72.0) go to 101
      if (x .lt.130.0 .or. x .gt.174.0) go to 101
      iaz = (iaz + roundr) / divsor
      x=x*(-1.)
      if(iaz.lt.1999)then
         ai=ai+1
         ax(ai)=x
         ay(ai)=y
         muv(ai)=iaz
      endif
      go to 101

  500 continue
C 
C******************************************************** 
C*    Plot selected lower 48 states stations            *
C********************************************************
C
      call mrclean(mark,ux,uy,iuv,ii)
      call mapstc ('OU - OUTLINE DATASET','US') 
      call gsplci (1)
      call maplot  
      call labtop (title, fecha, .013,0.0)
C 
C******************************************************** 
C*    Plot selected Hawaiian stations, if any           *
C********************************************************
C
      if (hi .gt. 0) then
         call gsclip(0)
         call gsasf(iasf)
         call mappos (0.,.30,.1,0.45 )  
         call maproj('ME',0.,-165.,0.)
         call mapset('CO - CORNERS',22.5,-161.,18.5,-154.) 
         call mapsti('GR',0)
         call mapsti('C7',0)
         call mapsti('C6',1)
         call mapstc ('OU - OUTLINE DATASET','CO') 
         call mapsti('C1 - PERIMETER',0)
         call mapint
         call gsplci (1)
c        sub country plots better outlines than ncar graphics
         call country
         ylatn = 22.5
         xlonw = -161.0
         ylats = 18.5
         xlone = -154.0
         xof   = .0005
         yof   = .001 
         call mrclean(mark,hx,hy,kuv,hi)
         call mapstc ('OU - OUTLINE DATASET','PS') 
         call gsplci (1)
         call labtop ('Hawaii', ' ', .03,.15)
      endif
C 
C************************************************* 
C*    Plot Puerto Rican stations, if any         *
C*************************************************
C
      if (pi .gt. 0) then
         call gsclip(0)
         call gsasf(iasf)
         call mappos (.50,.80,.13,0.33)  
         call maproj('ME',0.,-66.,0.)
         call mapset('CO - CORNERS',19.0,-67.5,17.5,-65.0) 
         call mapsti('GR',0)
         call mapsti('C7',0)
         call mapsti('C6',1)
         call mapstc ('OU - OUTLINE DATASET','CO') 
         call mapsti('C1 - PERIMETER',0)
         call mapint
         call gsplci (1)
c        sub country plots better outlines than ncar graphics
         call country
         ylatn = 19.0
         xlonw = -67.5
         ylats = 17.5
         xlone = -65.0
         xof   = -.0003
         yof   = .0005
         call mrclean(mark,px,py,luv,pi)
         call mapstc ('OU - OUTLINE DATASET','PS') 
         call gsplci (1)
         call labtop ('Puerto Rico', ' ', .02,.16)
      endif
C 
C********************************************** 
C*   Plot Alaskan stations, if any            *
C**********************************************
C
      if (ai .gt. 0) then
         call gsclip(0)
         call gsasf(iasf)
         call mappos (.63,.79,.70,0.85)  
c        call mappos (.82,.96,.31,0.46)  
         call maproj('ME',0.,-150.,0.)
         call mapset('CO',72.,-174.,54.,-130.)
         call mapsti('GR',0)
         call mapsti('C7',0)
         call mapsti('C6',1)
         call mapstc ('OU - OUTLINE DATASET','PO') 
         call mapsti('C1 - PERIMETER',0)
         call mapsti ('EL - ELLIPTICAL BOUNDARY',0)
         call mapint
         call gsplci (1)
         ylatn = 72.0
         xlonw = -174.0
         ylats = 54.0
         xlone = -130.0
         xof   = .01
         yof   = .01
         call mrclean(mark,ax,ay,muv,ai)
         call mapstc ('OU - OUTLINE DATASET','PO') 
         call gsplci (1)
         call maplot
         call labtop ('Alaska', ' ', .05,.14)
      endif
      call clsgks
      stop   
      end

      function dis(a,aa,b,bb)
      x=sqrt(((a-aa)**2)+((b-bb)**2) )
      dis=x
      return
      end

      subroutine mrclean(mark,zlon,zlat,izval,n)
      common/limits/ylatn,xlonw,ylats,xlone,xof,yof
      dimension zlon(n),zlat(n),izval(n)
      dimension idat(150,150),zdat(150,150)

      character *4 wet
      tdx=(xlone-xlonw)
      tdy=(ylatn-ylats)
      dx=tdx/150.
      dy=tdy/150.

c     init work arrays

      do i=1,150
         do j=1,150
            zdat(j,i)=-999
            idat(j,i)=-999
         enddo
      enddo

      do i=1,n
         lx=((zlon(i)-(xlonw))/dx)+1
         ly=((zlat(i)-(ylats))/dy)+1
c        write(*,477) lx,ly,zlon(i),zlat(i),izval(i)
 477     format (2i5,2f8.2,i5)
         if(zdat(lx,ly).lt.float(izval(i)))then
            idat(lx,ly)=i
            zdat(lx,ly)=float(izval(i))
            if(i.lt.10)then
c              write(*,666)i,lx,ly,zlat(i),zlon(i),izval(i)
  666          format(1x,3i5,2f8.1,i5)
            endif
         endif
      enddo

      dtest=dy
      do i=1,149
         do j=1,150
            m=idat(j,i)
            mm=idat(j,i+1)
            if(m.gt.0.and.mm.gt.0)then
               if((zlat(mm)-zlat(m)).lt.dtest)then
               endif
            endif
         enddo
      enddo

      do i=1,150
         do j=1,150
            if(idat(j,i).gt.0)then
               m=idat(j,i)
c              write(*,668)m,zlat(m),zlon(m),izval(m)
  668          format(1x,i5,2f8.2,i7)
               call maptrn(zlat(m),zlon(m),xloc,yloc)
               izzy=(izval(m))
               write(unit=wet,fmt='(i4)')izzy
               if(izval(m).lt.1999)then
                  if(izzy.gt.99)zilch=0.
                  if(izzy.lt.10)zilch=0.00234
                  if(izzy.ge.10.and.izzy.le.99)zilch=0.00117
                  if(izzy.lt.0)zilch=0.
                  zilch=zilch*0.6
                  call plchlq(xloc-zilch,yloc,wet,0.0035*0.6,0.,0.)
c                 write(*,344) xloc, yloc, zilch, wet
 344              format(3f10.4,a6)
                  if(izzy.ge.mark) then
c                    call plchlq(xloc+.075,yloc+.12,'o',.015*0.9,0.,0.)
                     call plchlq(xloc+xof,yloc+yof,'o',.015*0.9,0.,0.)
                  endif                             
               endif
            endif
         enddo
      enddo

      return
      end

      subroutine labtop (labl,labdat,size,yoff)
      character*(*) labl, labdat 

C Put a label just above the top of the plot.  The SET call is re-done
C to allow for the use of fractional coordinates, and the text extent
C capabilities of the package PLOTCHAR are used to determine the label
C position.
 
      call getset (xvpl,xvpr,yvpb,yvpt,xwdl,xwdr,ywdb,ywdt,lnlg)
      szfs=size*(xvpr-xvpl)
      call set (0.,1.,0.,1.,0.,1.,0.,1.,1)
      call pcgeti ('QU - QUALITY FLAG',IQUA)
      call pcseti ('QU - QUALITY FLAG',0) 
      call pcseti('CD',1)
      call pcseti ('TE - TEXT EXTENT COMPUTATION FLAG',1)
      call plchhq (.5,.5,labl,szfs,360.,0.)
      call pcgetr ('DB - DISTANCE TO BOTTOM OF STRING',DBOS)
      xloc=.45*(xvpl+xvpr)+0.01 
      szfs=szfs/1.5
      yloc= yvpt+szfs+DBOS 
      yyloc=yloc+(4.0*szfs)
      xxloc=xloc + (3.*szfs) 
      xszfs=szfs*6.0/5.0
      call plchhq (xloc,yyloc-0.01-yoff,labl,xszfs,0.,0.)  
      call plchhq (xloc,yloc-yoff,labdat,szfs*0.9,0.,0.)
                                                      
      B=0.4 *0.9
      call sflush 

      call pcseti ('QU - QUALITY FLAG',IQUA)
      call set    (xvpl,xvpr,yvpb,yvpt,xwdl,xwdr,ywdb,ywdt,lnlg)
      return 
      end

      subroutine country
      dimension x(5000),y(5000)
      character * 80 cstlin 
C  
C     this routine reads a file with better outlines for the 
C     islands of Hawaii and Puerto Rico than are in the NCAR system
C
      cstlin='/disk02/operations/cobprod/geography/coast-hi-pr.fil'
      irk=3
      call gsplci(1)
      call gstxci(1)
      call gslwsc(1.)
      call sflush
      open(10,file=cstlin,form='formatted',status='old')
 1    read(10,*,end=500,err=1)id
      i=1
 2    read(10,*,end=500,err=600)x(i),y(i)
      i=i+1
      go to 2
 600  i=i-1
      if(i.gt.1)then
 200  call mapit(y(1),x(1),0)
      do j=2,i
      call mapit(y(j),x(j),1)
      enddo
 3456 call mapiq
      go to 1
      endif
      go to 1
 500  irk=7
      call gslwsc(1.)
      call sflush
      close (10)
      return
      end