subroutine faavg(rlonf,rlatf,rlong,rlatg,rf1,rf2,rg,tg,
     +                 mr,ipert,tgavg,ierr)
c
c     This routine is for averaging a function of radius tg(rg)
c     centered at (rlong,rlatg) over a an annulus centered
c     at (rlonf,rlatf). The annulus is defined by rf1,rf2. 
c     
c     Input: 
c       rlonf,rlatf: Lat,lon in deg N and deg E of annulus center
c       rlong,rlatg: Lat,lon in deg N and deg E of center of the
c                    the function to be averaged.
c       rf1: Inner radius of annulus (km)
c       rf2: Outer radius of annulus (km)
c       rg(0:nr) - Array of radii (km) 
c       tgi(0:nr)- Array of function values to be averaged.
c                  Function values outside of rg(nr) are assumed
c                  to be zero. 
c       mr       - Dimension of rg and tgi
c       ipert    - Index of last radius to use for perturbation
c                  calculation. This is also the last radii
c                  that will be used in the area average. 
c                  If ipert < 1 this option is skipped.  
c     Output: 
c       tgavg    - annular average of tg
c       ierr     - Error flag, =0 for normal return, 
c                              =1 for error.
c
      dimension rg(0:mr),tg(0:mr)
c
c     Initialize values
      ierr  = 1
      tgavg = 0.0
      pi    = 3.14159265
      dtr   = pi/180.0
c     
c     Check input values
      if (rf1 .lt. 0.0) return
      if (rf2 .le. 0.0) return
      if (rf1 .ge. rf2) return
c
c     Calculate distance from average annulus center
c     to the function center.
      call distk(rlonf,rlatf,rlong,rlatg,dx,dy,d)
c
c     Calculate areas inside inner and outer radii of the annulus
      if (rf1 .gt. 0.0) then
         a1 = pi*rf1*rf1
      else
         a1 = 0.0
      endif
c
      a2 = pi*rf2*rf2
c
c     Determine outer radius of tg to include in the average,
c     and the constant for the perturbation calculation. 
      if (ipert .lt. 1) then
         nr = mr
         tgbar = 0.0
      else
         nr = ipert
         if (nr .gt. mr) nr = mr
         tgbar = tg(nr)
      endif
c
c     Calculate tgavg for the inner circular area
      tg1 = 0.0
      if (rf1 .gt. 0.0) then
         rf = rf1
         do i=1,nr
            ri  = rg(i)
            rim = rg(i-1)
c
            if (d .ge. rf) then
               if (abs(d-ri) .ge. rf) then
                  thetai  = 0.0
               else
                  if (ri .le. 0.0) then
                     thetai = 0.0
                  else
                     thetai = acos( (ri**2 + d**2 - rf**2)/(2.0*ri*d) )
                  endif
               endif
c
               if (abs(d-rim) .ge. rf) then
                  thetaim  = 0.0
               else
                  if (rim .le. 0.0) then
                     thetaim = 0.0
                  else
                     thetaim = acos( (rim**2 + d**2 - rf**2)/
     +                                          (2.0*rim*d) )
                  endif
               endif
c
            else
               if     (ri .gt. rf+d) then
                  thetai = 0.0
               elseif (ri .le. rf-d) then
                  thetai = pi
               else
                  if (ri .le. 0.0) then
                     thetai = 0.0
                  else
                     thetai = acos( (ri**2 + d**2 - rf**2)/(2.0*ri*d) )
                  endif
               endif
c
               if     (rim .ge. rf+d) then
                  thetaim = 0.0
               elseif (rim .le. rf-d) then
                  thetaim = pi
               else
                  if (rim .le. 0.0) then
                     thetaim = 0.0
                  else
                     thetaim = acos( (rim**2 + d**2 - rf**2)/
     +                                          (2.0*rim*d) )
                  endif
               endif
            endif
c
            tgmh = 0.5*(tg(i)+tg(i-1)) - tgbar 
            amh  = (ri**2-rim**2)*0.5*(thetai+thetaim)
c
            tg1 = tg1 + amh*tgmh
c
c        write(6,801) i,ri,rim,tg(i),tg(i-1),thetai/dtr,thetaim/dtr
c 801    format(i3,' ri=',f5.0,' ri-1=',f5.0,' ti=',f5.1,' ti-1=',f5.1,
c    +          ' thi=',f5.1,' thi-1=',f5.1) 
c
         enddo
         tg1 = tg1/a1
c
      endif
c
c     Calculate tgavg for the outer circular area
      tg2 = 0.0
      rf     = rf2
      do i=1,nr
         ri  = rg(i)
         rim = rg(i-1)
c
         if (d .ge. rf) then
            if (abs(d-ri) .ge. rf) then
               thetai  = 0.0
            else
               if (ri .le. 0.0) then
                  thetai = 0.0
               else
                  thetai = acos( (ri**2 + d**2 - rf**2)/(2.0*ri*d) )
               endif
            endif
c
            if (abs(d-rim) .ge. rf) then
               thetaim  = 0.0
            else
               if (rim .le. 0.0) then
                  thetaim = 0.0
               else
                  thetaim = acos( (rim**2 + d**2 - rf**2)/(2.0*rim*d) )
               endif
            endif
c
         else
            if     (ri .gt. rf+d) then
               thetai = 0.0
            elseif (ri .le. rf-d) then
               thetai = pi
            else
               if (ri .le. 0.0) then
                  thetai = 0.0
               else
                  thetai = acos( (ri**2 + d**2 - rf**2)/(2.0*ri*d) )
               endif
            endif
c
            if     (rim .ge. rf+d) then
               thetaim = 0.0
            elseif (rim .le. rf-d) then
               thetaim = pi
            else
               if (rim .le. 0.0) then
                  thetaim = 0.0
               else
                  thetaim = acos( (rim**2 + d**2 - rf**2)/(2.0*rim*d) )
               endif
            endif
         endif
c
         tgmh = 0.5*(tg(i)+tg(i-1)) - tgbar 
         amh  = (ri**2-rim**2)*0.5*(thetai+thetaim)
c
         tg2 = tg2 + amh*tgmh
c
c        write(6,802) i,ri,rim,tg(i),tg(i-1),thetai/dtr,thetaim/dtr
c 802    format(i3,' ri=',f5.0,' ri-1=',f5.0,' ti=',f5.1,' ti-1=',f5.1,
c    +          ' thi=',f5.1,' thi-1=',f5.1) 
c
      enddo
      tg2 = tg2/a2
c
c     Calculate the average tg over the annulus
      tgavg = (a2*tg2 - a1*tg1)/(a2-a1)
c
      ierr = 0
c
      return
      end