cdeck homer6
      subroutine homer6 (coor,n,npoly,area)
c   coor - 2 dimension array of coordinates
c   n    - number of coordinates in coor(n .ge. 3)
c   npoly - number of points in polygon, sorted in clockwise
c           order to top of coor on return (0 means error)
c   area  - area in polygon
       real coor(3,n), area
       integer n, npoly
      logical   posx
      include 'homerc.for'
c   zero   - smallest value allowed in denominator when checking slope
c   bigval - value of slope when denominator is zero
      real zero, bigval
      data zero /1.e-10/, bigval /1.e30/
c   find maximum y value
      ymax=coor(2,1)
      iymax=1
      do 10 i=2,n
         if (ymax.ge.coor(2,i)) go to 10
            ymax=coor(2,i)
            iymax=i
   10    continue
c   place maximum y coordinate in position 1
      call homer7 (coor,1,iymax)
      nm1=n-1
      np1=n+1
      do 80 i=1,nm1
        icurt=0
        slcurt=0.
        posx=.false.
        if (i.eq.1) alpha=0.
        if (i.ne.1) alpha=atan2(coor(2,i)-coor(2,i-1),coor(1,i)-coor(1,i
     1  -1))
         cosalp=cos(alpha)
         sinalp=sin(alpha)
         ip1=i+1
         do 70 k=ip1,np1
            if (k.ne.np1) go to 20
               if (i.eq.1) go to 70
               j=1
            go to 30
   20       j=k
c   coordinate axis rotation and translation such that coor(.,i)
c   is at the origin, with all other points in quadrats 2 and 3
c   (i.e., all other points 'below' the point in coor(.,i))
c   xp and yp are the transformed coordinates for coor(.,j)
   30       xp=(coor(1,j)-coor(1,i))*cosalp+(coor(2,j)-coor(2,i))*sinalp
            yp=(coor(2,j)-coor(2,i))*cosalp-(coor(1,j)-coor(1,i))*sinalp
            if (.not.posx) go to 40
               if (xp.lt.zero) go to 70
                  if ((yp/xp).gt.slcurt) go to 60
                  go to 70
   40       if (xp.lt.-zero) go to 50
               posx=.true.
               if (xp.gt.zero) go to 60
                  icurt=j
                  slcurt=bigval
                  go to 70
   50          if ((yp/xp).lt.slcurt) go to 70
   60             slcurt=yp/xp
                  icurt=j
   70    continue
         if (icurt.eq.1) go to 90
         call homer7 (coor,i+1,icurt)
   80 continue
   90 npoly=i
      area=0.
      if (npoly.lt.3) go to 110
         area=coor(1,1)*(coor(2,npoly)-coor(2,2))+coor(1,npoly)*(coor(2
     1   ,npoly-1)-coor(2,1))
         npolym1=npoly-1
         do 100 i=2,npolym1
            area=area+coor(1,i)*(coor(2,i-1)-coor(2,i+1))
  100       continue
         area=area*0.5
  110 return
      end