c    @(#)areai.f	1.2
c
c
c
      subroutine areai(a,iq1,iq2,iq3,a1,a2,a3)
c
c --- on output, a1 is the area average of array 'a'; a2 is the area
c --- average of the 'quad' operation on 'a'; a3 is from 'grd2' (short for
c --- 'grad squared').
c
      include 'qg_par.inc'
      common/cnn/n,nm1,nm2,np1
      common/parm/kz,m,mm1,mm2,mp1,msq,tim,tmax,tstart,hy,rea,dt,
     +  xbasin,yy(mx0),icalc,kzp1,kzm1,iftop,ifbot,iftvv,ifrel,theta,
     +   ifsbl,ifeva

      dimension a(mn0)
c
c --- function for the integration of (norm(del a))**2
c
      del(jp)=( (a(jp)-a(jp-1))**2 + 
     +          (a(jp)-a(jp-1))*(a(jp-m)-a(jp-m-1)) +
     +          (a(jp-m)-a(jp-m-1))**2 +
     +          (a(jp)-a(jp-m))**2 +
     +          (a(jp)-a(jp-m))*(a(jp-1)-a(jp-m-1)) +
     +          (a(jp-1)-a(jp-m-1))**2)/3
c
      a1=0.
      a2=0.
      a3=0.
c
      ml=mm1/4 + 1
      mlp1=ml+1
      mu=m-ml + 1
      mln=nm1/4+1
      mlp1n=mln+1
      mun=n-mln+1
c
c --- first integrate a and a squared.
c
      do 5 j=mln,mun
        ip=(j-1)*m  +  ml-1
        do 4 i=ml,mu
          ip=ip+1
          if(iq1.ne.0) a1 = a1 + a(ip)
          if(iq2.ne.0) a2 = a2 + a(ip)*a(ip)
    4   continue
    5 continue
c
c --- now do norm of del squared.
c
      do 20 j=mlp1n,mun
        ip=(j-1)*m  +  ml
        do 10 i=mlp1,mu
          ip=ip+1
          if(iq3.ne.0) a3=a3+del(ip)
   10   continue
   20 continue
c
c --- normalize
c
      nmul=mu-ml+1
	nmuln=mun-mln+1
	xnm=nmuln*nmul
      a1=a1/xnm
      a2=a2/xnm
      a3=a3/(hy*hy*xnm)
      return
      end
