cdeck homr3a
      subroutine homr3a (xtfm,ytfm,qa,ndim,storxy,nintmax,temp)
      integer nintmax,ndim
      real qa(nintmax,nintmax), xtfm(nintmax), ytfm(nintmax), storxy(3
     1 ,ndim), temp(1)
      include 'maincomo.for'
      include 'homerc.for'
c     qa is the two dimensional frequency array
c     xtfm is the array of upper bounds for the x axis
c     ytfm is the array of upper bounds for the y axis
      character *10 answer(4)
      character *3 yesno(2)
      integer fdate(5)
      logical   qfinf, qfing
      data answer /'no map','fixes only','fill also','both ways'/
      data yesno /'no','yes'/
c
c   read first case in the subfile
c
   10 nrt=0
      call homer2 (pxfix,pyfix,ptime,qfing,qfinf,temp)
      if (qfinf) return
      if (qfing) go to 400
      nrt=1
c     zero out all arrays and fill xtfm and ytfm
c   set flag for first fix outside of map rectangle
      nerst=0
      nherb=0
      nsqt=0
      nsqtv=0
      nkim=ndim
c     with nomap=1, value of ysiz is slightly larger than old version
      ysiz=siz(idic)+siz(idic)/9.0*nomap
      naddf=0
      nzim=1
      if (nint(idic).lt.1) nint(idic)=100
      do 20 i=1,nint(idic)
      do 20 j=1,nint(idic)
   20 qa(i,j)=0
      ierror=0
      xmax=-0.1e-30
      xmin=+0.1e+30
      ymax=-0.1e-30
      ymin=+0.1e+30
      xtfm(1)=xco(idic)+siz(idic)
      ytfm(1)=yco(idic)+ysiz
      do 30 i=2,nint(idic)
      xtfm(i)=xtfm(i-1)+siz(idic)
   30 ytfm(i)=ytfm(i-1)+ysiz
c     set interval variables for use at statement 18
      xxxco=xco(idic)-.00001
      yyyco=yco(idic)-.00001
   40 if (pxfix.gt.xco(idic).and.pxfix.lt.xtfm(nint(idic)).and.pyfix.gt
     1 .yco(idic).and.pyfix.lt.ytfm(nint(idic))) go to 50
      write (iolp,430) pxfix,pyfix
      call linect (1)
      call homer2 (pxfix,pyfix,ptime,qfing,qfinf,temp)
      if (qfing) go to 400
      go to 40
   50 nzim=0
      do 60 i=1,nint(idic)
      if (pxfix.lt.xtfm(i)) go to 70
   60 continue
   70 ip=i
      do 80 i=1,nint(idic)
      if (pyfix.lt.ytfm(i)) go to 90
   80 continue
   90 jp=i
      qa(ip,jp)=1
      nallf=1
c
c   set date and time of first case used
c
      call homer9 (fdate,ptime)
      if (nkim.eq.0) go to 100
      storxy(1,1)=pxfix
      storxy(2,1)=pyfix
      storxy(3,1)=ptime
c
c     loop 140 contains vector fill, placement of fixes in qa matrix,
c     and storage of fixes for homer5
c
  100 call homer2 (x,y,atime,qfing,qfinf,temp)
      if (qfing) go to 150
      nrt=nrt+1
c     are the coordinates within the proper range
      if (x.lt.xco(idic).or.xtfm(nint(idic)).lt.x.or.y.lt.yco(idic).or
     1 .ytfm(nint(idic)).lt.y) then
         xmin=min(x,xmin)
         ymin=min(y,ymin)
         xmax=max(x,xmax)
         ymax=max(y,ymax)
         ierror=ierror+1
         if (nerst.eq.0) then
            nerst=nallf
            ndon(idic)=0
         endif
      endif
      ii=(x-xxxco)/siz(idic)+1.0
      jj=(y-yyyco)/ysiz+1.0
c     above is faster method of determining ii and jj than searching
c     arrays xtfm and ytfm - need to check roundoff correction  ********
c
c   vector fill section
c
c   if nerst not set to zero, then fix out of bounds has occured
c   no more vector fill after the nerst fix
      if (nerst.ne.0) go to 140
c   check to see if vectime parameter is exceeded
      if (atime.ge.(ptime+time(idic))) go to 140
c
c   check to see if vecdist parameter exceeded
c
      if (sqrt((x-pxfix)**2+(y-pyfix)**2).gt.tobig(idic)) go to 140
c
c   no fill needed if both locations in same cell
c
      if (ii.eq.ip.and.jj.eq.jp) then
c
c   horizonal vector fill
c
      elseif (ii.ne.ip.and.jj.eq.jp) then
         nst=min(ii,ip)
         nfin=max(ii,ip)
         ndiff=nfin-nst-2
         if (ndiff+1.gt.0) then
            nst=nst+1
            nfin=nst+ndiff
            do 110 kk=nst,nfin
            naddf=naddf+1
  110       qa(kk,jj)=qa(kk,jj)+1.0e-6
         endif
c
c   vertical vector fill
c
      elseif (ii.eq.ip.and.jj.ne.jp) then
         nst=min(jj,jp)
         nfin=max(jj,jp)
         ndiff=nfin-nst-2
         if (ndiff+1.gt.0) then
            nst=nst+1
            nfin=nst+ndiff
            do 120 kk=nst,nfin
            naddf=naddf+1
  120       qa(ii,kk)=qa(ii,kk)+1.0e-6
         endif
c
c     both the i and j componts are not equal -- have to calc a vector
c
      else
         nst=iabs(ip-ii)*iabs(jp-jj)
         if (nst.eq.1) go to 140
         xdiff=x-pxfix
         ydiff=y-pyfix
         ang=atan2(ydiff,xdiff)
         qxfix=pxfix
         qyfix=pyfix
         amag=sqrt(xdiff*xdiff+ydiff*ydiff)
         xl=siz(idic)*cos(ang)
         yl=siz(idic)*sin(ang)
  130    xnew=pxfix+xl
         ynew=pyfix+yl
         xot=xnew-qxfix
         yot=ynew-qyfix
         bmag=sqrt(xot*xot+yot*yot)
         if (bmag.ge.amag) go to 140
         inew=(xnew-xxxco)/siz(idic)+1.0
         jnew=(ynew-yyyco)/ysiz+1.0
c        above is faster method of determining indices than searching
c        arrays xtfm and ytfm - need to check roundoff correction  *****
         if (jnew.ne.jp.or.inew.ne.ip) then
            qa(inew,jnew)=qa(inew,jnew)+1.0e-6
            naddf=naddf+1
            ip=inew
            jp=jnew
         endif
         pxfix=xnew
         pyfix=ynew
         go to 130
      endif
c
c   increment number of fixes in cell ii,jj
c
  140 qa(ii,jj)=qa(ii,jj)+1
      pxfix=x
      pyfix=y
      ip=ii
      jp=jj
      ptime=atime
      nallf=nallf+1
      if (nkim.eq.0) go to 100
      if (nallf.gt.nkim) then
         nkim=0
         write (iolp,440)
      call linect (1)
      else
         storxy(1,nallf)=x
         storxy(2,nallf)=y
         storxy(3,nallf)=atime
      endif
      go to 100
c
c   a subfile of data has been read
c
  150 call dphead
      write (unit=iolp,fmt=410) (varnam(varnum(i)),(varlab(j,i),j=1,5),i
     1 =1,3)
      call linect (7)
      write (iolp,420) xco(idic),yco(idic),siz(idic),nint(idic),time
     1 (idic),tobig(idic),ndon(idic),answer(nsum+1),yesno(nomap+1),yesno
     2 (irem+1),qoptsr(3),qoptsr(4),qoptsr(6)
      call linect (13)
      if (qoptsr(7).or.qoptsr(8).or.qoptsr(9).or.qoptsr(10).or.qoptsr(11
     1 )) write (unit=iobcd2,fmt='(1h*,a8,1h*)') sfdic(idic)
      nzim=naddf
      if (nrt.lt.4) go to 400
      naddf=0
      if (ierror.eq.0) go to 160
      write (iolp,450) ierror,nerst,xmin,xmax,ymin,ymax
      call linect (5)
      write (iolp,460)
      call linect (1)
c
c     determine bounds of submatrix containing fixes
c
c   matrix inside rectange will be defined by indices
c   x -- ii-iend      y -- jj-jend
  160 do 170 i=1,nint(idic)
      do 170 j=1,nint(idic)
      if (qa(i,j).gt.0.) go to 190
  170 continue
  180 write (iolp,470)
      call linect (2)
      return
  190 ii=i-1
      if (ii.le.0) ii=1
      do 200 j=1,nint(idic)
      do 200 i=1,nint(idic)
      if (qa(i,j)) 210,200,210
  200 continue
      go to 180
  210 jj=j-1
      if (jj.le.0) jj=1
      do 220 i=1,nint(idic)
      ki=nint(idic)-i+1
      do 220 j=1,nint(idic)
      if (qa(ki,j)) 230,220,230
  220 continue
      go to 180
  230 iend=ki+1
      if (iend.gt.nint(idic)) iend=nint(idic)
      do 240 j=1,nint(idic)
      kj=nint(idic)-j+1
      do 240 i=1,nint(idic)
      if (qa(i,kj)) 250,240,250
  240 continue
      go to 180
  250 jend=kj+1
      if (jend.gt.nint(idic)) jend=nint(idic)
      if (ndon(idic).eq.0) go to 350
c
c     home range fill section
c
  260 nsmor=0
c     search for an initial non zero value in a column
      do 300 j=jj,jend
         do 290 k=ii,iend-ndon(idic)
            if (qa(k,j).gt.0.0) then
c              an initial non zero point found
               do 280 kk=k+ndon(idic),k+1,-1
c              search for non negative qa in the string
                  if (qa(kk,j).gt.0.) then
                     do 270 kkk=k+1,kk-1
                        if (qa(kkk,j).eq.0.0) then
                           qa(kkk,j)=+1.0e-6
                           nsmor=nsmor+1
                        endif
  270                   continue
                     go to 290
                  endif
  280             continue
               go to 290
            endif
  290       continue
  300    continue
c
c   search for an initial non zero value in a row
c
      do 340 i=ii,iend
         do 330 k=jj,jend-ndon(idic)
            if (qa(i,k).gt.0.0) then
c              an initial non zero point found
               do 320 kk=k+ndon(idic),k+1,-1
c               search for non negative qa in the string
                  if (qa(i,kk).gt.0.) then
                     do 310 kkk=k+1,kk-1
                        if (qa(i,kkk).eq.0.0) then
                           qa(i,kkk)=+1.0e-6
                           nsmor=nsmor+1
                        endif
  310                   continue
                     go to 330
                  endif
  320             continue
               go to 330
            endif
  330       continue
  340    continue
      naddf=naddf+nsmor
c   the fill option is recursive, so take another shot
c   if some squares were filled this time
      if (nsmor.gt.0) go to 260
  350 do 360 i=1,nint(idic)
      do 360 j=1,nint(idic)
         if (qa(i,j).ge.1.0) then
            nherb=nherb+1
         else
            if (qa(i,j).gt.0.0) nsqt=nsqt+1
         endif
  360    continue
      nsqtv=nsqt-naddf
      notin=0
      if (irem.eq.0) go to 390
c     wild fix removal section - note wild fixes are presently included
c     minimum area method calculation even if deleted here. to remedy th
c     must purge storxy of wild fixes, shorten arrays appropriatly
c     and adjust nallf
      do 380 i=1,nint(idic)
      do 380 j=1,nint(idic)
      if (qa(i,j).eq.0) go to 380
      ist=max(i-1,1)
      jst=max(j-1,1)
      ibox=min(i+1,nint(idic))
      jbox=min(j+1,nint(idic))
      nber=0
      do 370 ir=ist,ibox
      do 370 jr=jst,jbox
      if (qa(ir,jr).gt.0) nber=nber+1
      if (nber.gt.1) go to 380
  370 continue
      qa(i,j)=0
      notin=notin+1
  380 continue
  390 ltot=nallf+naddf+nzim
      ktot=nherb+nsqt-notin
      sqmeter=ktot*siz(idic)*ysiz
      hect=sqmeter*0.0001
      acres=sqmeter*0.0002471054
      write (iolp,480) fdate
      call linect (3)
      call homer9 (fdate,ptime)
      write (iolp,490) fdate
      call linect (3)
      write (iolp,500) nrt,nallf,nherb,naddf,nzim,nsqtv,nsqt,notin,ltot
      call linect (10)
      if (notin.gt.0) then
         write (iolp,510)
         call linect (2)
      endif
      write (iolp,520) ktot,hect,acres
      call linect (2)
c
c   perform minimum area polygon, dunn, and jennrich-turner
c   home range estimates
c
      if (nkim.gt.0) then
         call homer5 (nallf,ndim,storxy)
      endif
      if (nsum.gt.0) then
c
c   print home range maps
c
         call homer4 (qa,xtfm,ytfm,ii,jj,iend,jend,nintmax)
      endif
      go to 10
c
c   error return for too few cases in a subfile
c
  400 write (unit=iolp,fmt=530) nrt
      call linect (3)
      go to 10
c
  410 format (///10x,'v a r i a b l e   s u m m a r y'/10x,'x coordinate
     1 - ',a8,1h,,5a8/10x,'y coordinate - ',a8,1h,,5a8/10x,'time
     2 - ',a8,1h,,5a8)
  420 format (/10x,'value in meters of x lower bound (xorigin=) is ',g14
     1 .8/10x,'value in meters of y lower bound (yorigin=) is ',g14.8/10
     2 x,'size of an interval (square=) is ',g14.8,' meters'/10x,
     3'number of intervals (grids=) is ',i4/10x,
     4 'maximum time (vectime=) between fixes for fill ',g14.8,
     5' days'/10x,'distance limit for vector fill (vecdist=) is ',
     1   g14.8,' meters'
     6   /10x,'home range fill (fill=) option is ',i4,' squares'
     7   /10x,'home range map print option (map=) is ',a
     8   /10x,'map overlay option (mapscale) is ',a
     9   /10x,'wild fix removal option (remove) is ',a
     a   /10x,'no minimum area method option (option 3) is ',l4
     b   /10x,'no jennrich-turner home range option (option 4) is ',l4
     c   /10x,'dunn home range option (option 6) is ',l4)
  430 format (10x,'fix outside range specified - it has been discarded
     1  x = ',g14.8,3x,'y = ',g14.8)
  440 format (10x,'no minimum area method calculations for this run - in
     1crease value of maxcase to no. of fixes in subfile')
  450 format (/10x,'there were ',i6,' cases with fixes outside the range
     1 specified.'/10x,' there were were no fixes added after the ',i6,'
     2st case.,/10x,set your interval to include  these values.  ignore
     3 a value of +0.1e+30',/10x,'xmin=',g14.8,5x,'xmax=',g14.8,5x,'ymin
     4=',g14.8,5x,'ymax=',g14.8)
  460 format (10x,'home range table printed does not include  fixes added
     1 ')
  470 format (10x,'error, no fixes in specified range - check data input
     1 and control card parameters'/10x,'no home range output for this d
     2ata set')
  480 format (/10x,'date and time on first case used - ',2(i2.2,1h/),i2.
     1 2,1x,2i2.2)
  490 format (/10x,'date and time on last case used  - ',2(i2.2,1h/),i2.
     1 2,1x,2i2.2)
  500 format (/10x,'number of fixes in subfile',34x,i5/10x,'number of fi
     1xes (from subfile) used',25x,i5/10x,'number of squares filled by a
     2nimal fixes',20x,i5/10x,'number of fixes added = squares filled by
     3 home range fill',3x,i5/10x,'number of fixes added by tracking-vec
     4tor fill',15x,i5/10x,'number of squares added by tracking-vector f
     5ill',13x,i5/10x,'total squares added by fill options',25x,i5/10x,'
     6number of squares removed by wild data option',15x,i5/10x,'total n
     7umber of fixes in home range printout',16x,i5)
  510 format (10x,'note - numbers of fixes actually used after removal b
     1y wild data option unknown'/17x,'minimum area method calculations
     2 include  fixes in squares removed by wild data option')
  520 format (10x,'total number of squares in home range printout',14x
     1 ,i5/10x,'grid cell home range  ',g14.8,' hectares   (',g14.8,' ac
     2res).')
  530 format (/10x,'* * *   e r r o r   * * *   only',i2,' cases in subf
     1ile.'/10x,'processing of this subfile is skipped.')
      end