c Extracts plouff stations from INSIDE. a polygonal geographic border.
c Also saves outside ones in a separate file "OUTSIDE.DAT"
      character a*1,b*1,sta*8,rest*57
      common /body/ xb(201),yb(201),nb
      common /coord/ latcen,loncen,cenlat,cenlon
      data nc,ni,no,iw,n,ifdel/6*0/,x0,y0/2*0.0/
      call prom (ifdel)
      if (ifdel .eq. 999) stop
      if (ifdel .eq. 99) go to 999
      n1=nb-1
    1 n=n+1
      read (15,500,err=9,end=99) sta,a,lat,jlat,b,lon,jlon,rest
  500 format (bz,a8,a1,i2,i4,a1,i3,i4,a57)
      dp=(lat-latcen)+(0.01*jlat-cenlat)/60.
      dl=(loncen-lon)+(cenlon-0.01*jlon)/60.
      call proj (dp,dl,x0,y0)
      call wind (iw,n1,x0,y0)
      if (iw) 2,3,4
c Outside polygon
    2 no=no+1
      if (ifdel .ne. 0) go to 1
      write (9,500) sta,a,lat,jlat,b,lon,jlon,rest
      go to 1
c On border
    3 flat=0.01*jlat
      flon=0.01*jlon
      write (16,600) n,sta,lat,flat,lon,flon
  600 format (' line',i5,', station',a8,i3,f6.2,i4,f6.2,' on border')
      nc=nc+1
c Inside polygon
    4 write (7,500) sta,a,lat,jlat,b,lon,jlon,rest
      ni=ni+1
      go to 1
    9 print 601, n
      write (16,601) n
  601 format (' Stop. Format error on line',i5,/,' Your output',
     1 ' file will terminate before it.')
   99 n=n-1
      print 609, n,ni,nc
      write (16,609) n,ni,nc
  609 format (i6,' input stations',/,i6,' stations inside or on ',
     1 'boundary are in file INSIDE.DAT',/,' (',i6,' stations on ',
     2 'border)',/,' You may want to delete trailing blank spaces ',
     3 'from INSIDE.DAT.')
      if (ifdel .eq. 1) go to 999
      close (9)
      print 610, no
      write (16,610) no
  610 format (i6,' stations outside the border are in a file named ',
     1 'OUTSIDE.DAT.')
  999 print 900
  900 format (' Print file SURROUND.PNT for a record of this session.')
      close (15)
      close (16)
      close (7)
      close (8)
      stop
      end
      SUBROUTINE PROJ(DPI,DLI,DP,DL)
      COMMON /PARAM/ AGE,AGEF,CP,SP,FM1,SC228,SC5
      DPR=1.745329E-2*DPI
      DLB=1.745329E-2*DLI
C OUTPUT Y-COORDINATE IS DP. X IS DL. BOTH ARE IN KILOMETERS.
C INPUT DL IS CENTRAL MERIDIAN OF MAP MINUS LONGITUDE IN DEGREES.
C INPUT DP IS LATITUDE MINUS CENTRAL LATITUDE OF MAP, IN DEGREES.
      G=(AGE+DPR*AGEF)*DLB
      DL=       G*(CP-DPR*SP)
      DP=        DPI*(FM1+DPR*SC228)+    DLB*G*SC5
      RETURN
      END
      SUBROUTINE PROJ3 (DPI)
      COMMON /PARAM/ AGE,AGEF,CP,SP,FM1,SC228,SC5
      DATA A,EQ,DEGR/6378.206, 6.768658E-3,1.745329E-2/
      DPR=DEGR*DPI
      CP=COS(DPR)
      SP=SIN(DPR)
      CP2=CP*CP
      SP2=1.0-CP2
      SCP=SP*CP
      SC5=0.5*SCP
      F1141=1.1417-9.6E-3*CP2
      FM1=(111.6993-CP2*F1141)*0.9996
      SC228=F1141*SCP*0.9996
      GE=1.0+0.5*EQ*SP2*(1.0+0.75*EQ*SP2)
      AGE=A*GE*0.9996
      AGEF=A*EQ*SCP*0.9996
      WRITE (16,200)
  200 FORMAT (10X,'UNIVERSAL TRANSVERSE MERCATOR PROJECTION')
      RETURN
      END
      subroutine wind(iw,n1,x0,y0)
c finds winding number. 4 if inside polygon. -1 outside.
c  zero on the boundary (iw=). x,y corners. x0,y0 point.
c Requires that last point is a duplicate of first; n1=total number
c of sides of the polygon (1 less than required points).
c Algorithm from Alan Cogbill based on Leland, K, 1975, Math of
c  Computation, v. 29, p.554-558.
      common /body/ x(201),y(201),n
      is=0
      x2=x(1)-x0
      y2=y(1)-y0
      k=0
   10 k=k+1
      if (k .gt. n1) go to 60
      ip=0
      it=0
      x1=x2
      y1=y2
      x2=x(k+1)-x0
      y2=y(k+1)-y0
      if (y2 .ge. 0.0) go to 20
c y2<0
      if (y1 .lt. 0.0) go to 30
c y1>=0, y2<0
      ip=1
      go to 40
   20 if (y1 .ge. 0.0) go to 30
c y2>=0, y1<0
      ip=-1
      go to 40
c end step1, begin step 2
   30 if (y1 .ne. 0.0) go to 10
      if (y2 .ne. 0.0) go to 10
      if (x1 .eq. 0.0) go to 80
      if (x2 .eq. 0.0) go to 80
      if (sign(1.0,x1) .ne. sign(1.0,x2)) go to 80
      go to 10
c end of step 2, begin step 3
   40 x1y2=x1*y2
      det=x2*y1-x1y2
c test value depends on computer
      if (abs(det) .le. 5.0e-7*abs(x1y2)) go to 80
      if (float(ip)*det .gt. 0.0) go to 50
      it=-ip
   50 is=is+it
      go to 10
   60 if (is .eq. 0) go to 70
      iw=iabs(is)
      return
   70 iw=-1
      return
   80 iw=0
      return
      end
      subroutine prom (ifdel)
      dimension latb(201),lonb(201),gtb(201),gnb(201)
      character infile*65,polfile*65,yes*1,clock*8,dat*9
      logical exists
      common /body/ xb(201),yb(201),nb
      common /coord/ latcen,loncen,cenlat,cenlon
      data tmax,emax,tmin,emin /2*0.0,90.,400./
      print 100
  100 format (' SURROUND, Plouff.',/,' To extract stations from ',
     1 'inside and outside a polygonal geographic boundary.')
      inquire (file='SURROUND.PNT',exist=exists)
      if (exists) then
        print 101
  101   format (' A previous print file SURROUND.PNT already exists.',
     1   /,' Do you want to STOP to save it?')
        read 501, yes
  501   format (a1)
        if (yes .eq. 'y' .or. yes .eq. 'Y' .or. yes .eq. ' ') then
          ifdel=999
          return
        end if
      end if
      inquire (file='INSIDE.DAT',exist=exists)
      if (exists) then
        print 102
  102   format (' A plouff output file INSIDE.DAT from this program ',
     1   'already exists.',/,' Do you want to STOP to save it?')
        read 501, yes
        if (yes .eq. 'y' .or. yes .eq. 'Y' .or. yes .eq. ' ')  then
          ifdel=999
          return
        end if
      end if
      print 103
  103 format (' TYPE the name of the plouff input file:')
      read 565, infile
  565 format (a65)
      inquire (file=infile,exist=exists)
      if (.not. exists) then
        print 104
  104   format (' File was not found.  Try once more.')
        print 103
        read 565, infile
        inquire (file=infile,exist=exists)
        if (.not. exists) then
          print 105
  105     format (' **STOP. File was not found.')
          ifdel=999
          return
        end if
      end if
      open (16,file='SURROUND.PNT',form='formatted',status='unknown')
      call date(dat)
      call time(clock)
      write (16,603) dat,clock,infile
  603 format (' SURROUND, Plouff, ',a9,1x,a8,/,' To extract stations ',
     1 'from inside or outside a polygonal geographic boundary.',/,
     2 ' Input file: ',a65)
      open (15,file=infile,form='formatted',status='old',blank='zero')
      print 106
  106 format (' You need a file with consecutive polygon corners.  ',
     1 'Each line has 4 decimal',/,'  numbers (separated by spaces or',
     2 ' commas) consisting of longitude (d,m)',/,'  followed by ',
     3 'latitude (d,m).  The first vertex need not be repeated.')
      print 107
  107 format (' TYPE the name of the polygon file:')
      read 565, polfile
      inquire (file=polfile,exist=exists)
      if (.not. exists) then
        print 104
        print 107
        read 565, polfile
        inquire (file=polfile,exist=exists)
        if (.not. exists) then
          print 108
          write (16,108)
  108     format (' **STOP. Polygon file was not found.')
          ifdel=999
          close (16)
          return
        end if
      end if
      write (16,205) polfile
  205 format (' Your file of polygon corners is:',/,3x,a65,/,' Your ',
     1 'file of stations inside the polygon is named INSIDE.DAT')
      ifdel=1
      nb=0
      print 203
  203 format (' Do you want the exterior stations to be saved in a ',
     1 'file?')
      read 501, yes
      if (yes .eq. 'y' .or. yes .eq. 'Y' .or. yes .eq. ' ') then
        ifdel=0
        print 201
        write (16,201)
  201   format (' The file for stations requested outside the polygon',
     1   ' is named OUTSIDE.DAT')
        inquire (file='OUTSIDE.DAT',exist=exists)
        if (exists) then
          print 109
  109     format (' A plouff output file OUTSIDE.DAT from this prog',
     1     'ram already exists.',/,' Do you want to STOP to save it?')
          read 501, yes
          if (yes .eq. 'y' .or. yes .eq. 'Y' .or. yes .eq. ' ')  then
            write (16,160)
  160       format (' You chose to STOP to save file OUTSIDE.DAT')
            ifdel=999
            close (16)
            return
          end if
        end if
        open (9,file='OUTSIDE.DAT',form='formatted',status='unknown')
      end if
      open (15,file=infile,form='formatted',status='old',blank='zero')
      open (8,file=polfile,form='formatted',status='old',blank='zero')
      open (7,file='INSIDE.DAT',form='formatted',status='unknown')
    4 read (8,*,end=7) gn,gm,ft,fm
c Change to positive coordinates consistent with plouff format
      gn=abs(gn)
      ft=abs(ft)
c Delimiter -999,0,-999,0 of usual other-lines file not needed
      if (gn .ge. 360 .or. ft .ge. 90) go to 7
      nb=nb+1
      if (nb .lt. 201) go to 9
      print 210
      write (16,210)
  210 format (' Stop. 200 vertices exceeded.')
   12 ifdel=99
      return
    9 flon=gn+gm/60.
      flat=ft+fm/60.
      if (flon .gt. emax)emax=flon
      if (flat .gt. tmax) tmax=flat
      if (flat .lt. tmin)tmin=flat
      if (flon .lt. emin) emin=flon
      latb(nb)=flat
      lonb(nb)=flon
      gtb(nb)=(flat-latb(nb))*60.0
      gnb(nb)=(flon-lonb(nb))*60.0
      write (16,602) nb,latb(nb),gtb(nb),lonb(nb),gnb(nb)
  602 format (' Corner:',i4,2(i4,'D',f6.2,'M'))
      go to 4
    7 print 208,nb
      write (16,208) nb
  208 format (i4,' corner coordinates read.')
      if (nb .gt. 2) go to 14
   13 print 202
      write (16,202)
  202 format (' ***STOP. Too few corners.')
      go to 12
   14 clat=0.5*(tmax+tmin)
      cen=0.5*(emax+emin)
      latcen=clat
      cenlat=60.*(clat-latcen)
      loncen=cen
      cenlon=60.*(cen-loncen)
      print 209, latcen,cenlat,loncen,cenlon
      write (16,209) latcen,cenlat,loncen,cenlon
  209 format (' Area center at',i3,'D',f5.1,'M latit',i4,'D',f5.1,
     1 'M longit')
      call proj3 (clat)
      do 8 j=1,nb
      dp=(latb(j)-latcen)+(gtb(j)-cenlat)/60.
      dl=(loncen-lonb(j))+(cenlon-gnb(j))/60.
    8 call proj(dp,dl,xb(j),yb(j))
c Winding algorithm requires last point is repeat of first point
      if (xb(1) .eq. xb(nb) .and. yb(1) .eq. yb(nb)) then
        if (nb .eq. 3) go to 13
        return
      end if
      nb=nb+1
      xb(nb)=xb(1)
      yb(nb)=yb(1)
      write (16,204)
  204 format (' A repeat of the first point was added to close the ',
     1 'polygon.')
      return
      end