ecearth3/sources/oasis3-mct/lib/mct/mct/m_SpatialIntegral.F90

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!-----------------------------------------------------------------------
! CVS m_SpatialIntegral.F90,v 1.22 2008-05-12 01:59:33 jacob Exp
! CVS MCT_2_8_0 
!BOP -------------------------------------------------------------------
!
! !MODULE: m_SpatialIntegral - Spatial Integrals and Averages using a GeneralGrid
!
! !DESCRIPTION:  This module provides spatial integration and averaging 
! services for the MCT.  For a field $\Phi$ sampled at a point ${\bf x}$ 
! in some multidimensional domain $\Omega$, the integral $I$ of 
! $\Phi({\bf x})$ is
! $$ I = \int_{\Omega} \Phi ({\bf x}) d\Omega .$$ 
! The spatial average $A$ of $\Phi({\bf x})$ over $\Omega$ is
! $$ A = {{ \int_{\Omega} \Phi ({\bf x}) d\Omega} \over 
! { \int_{\Omega} d\Omega} }. $$
! Since the {\tt AttrVect} represents a discretized field, the integrals 
! above are implemented as:
! $$ I = \sum_{i=1}^N \Phi_i \Delta \Omega_i $$
! and 
! $$ A = {{\sum_{i=1}^N \Phi_i \Delta \Omega_i } \over 
!{\sum_{i=1}^N \Delta \Omega_i } }, $$
! where $N$ is the number of physical locations, $\Phi_i$ is the value 
! of the field $\Phi$ at location $i$, and $\Delta \Omega_i$ is the spatial 
! weight (lenghth element, cross-sectional area element, volume element, 
! {\em et cetera}) at location $i$.
!
! MCT extends the concept of integrals and area/volume averages to include 
! {\em masked} integrals and averages.  MCT recognizes both {\em integer}
! and {\em real} masks.  An integer mask $M$ is a vector of integers (one
! corresponding to each physical location) with each element having value 
! either zero or one.  Integer masks are used to include/exclude data from
! averages or integrals.  For example, if one were to compute globally 
! averaged cloud amount over land (but not ocean nor sea-ice), one would 
! assign a $1$ to each location on the land and a $0$ to each non-land 
! location.  A {\em real} mask $F$ is a vector of real numbers (one corresponding 
! to each physical location) with each element having value within the 
! closed interval $[0,1]$.  .Real masks are used to represent fractional 
! area/volume coverage at a location by a given component model.  For 
! example, if one wishes to compute area averages over sea-ice, one must 
! include the ice fraction present at each point.  Masked Integrals and 
! averages are represented in the MCT by:
! $$ I = \sum_{i=1}^N {\prod_{j=1}^J M_i} {\prod_{k=1}^K F_i} 
! \Phi_i \Delta \Omega_i $$
! and 
! $$ A = {{\sum_{i=1}^N \bigg({\prod_{j=1}^J M_i}\bigg) \bigg( {\prod_{k=1}^K F_i}
! \bigg) \Phi_i 
! \Delta \Omega_i } \over 
!{\sum_{i=1}^N \bigg({\prod_{j=1}^J M_i}\bigg) \bigg( {\prod_{k=1}^K F_i} \bigg) 
!  \Delta \Omega_i } }, $$
! where $J$ is the number of integer masks and $K$ is the number of real masks.
!
! All of the routines in this module assume field data is stored in an 
! attribute vector ({\tt AttrVect}), and the integration/averaging is performed 
! only on the {\tt REAL} attributes.  Physical coordinate grid and mask 
! information is assumed to be stored as attributes in either a 
! {\tt GeneralGrid}, or pre-combined into a single integer mask and a single 
! real mask.  
!
! !INTERFACE:

 module m_SpatialIntegral
     
      implicit none

      private	! except

! !PUBLIC MEMBER FUNCTIONS:

      public :: SpatialIntegral        ! Spatial Integral
      public :: SpatialAverage         ! Spatial Area Average

      public :: MaskedSpatialIntegral  ! Masked Spatial Integral
      public :: MaskedSpatialAverage   ! MaskedSpatial Area Average

      public :: PairedSpatialIntegrals ! A Pair of Spatial 
                                      ! Integrals 

      public :: PairedSpatialAverages  ! A Pair of Spatial 
                                      ! Area Averages

      public :: PairedMaskedSpatialIntegrals ! A Pair of Masked
                                             ! Spatial Integrals 

      public :: PairedMaskedSpatialAverages ! A Pair of Masked
                                            ! Spatial Area Averages

      interface SpatialIntegral ; module procedure &
	   SpatialIntegralRAttrGG_
      end interface
      interface SpatialAverage ; module procedure &
	   SpatialAverageRAttrGG_ 
      end interface
      interface MaskedSpatialIntegral ; module procedure &
	   MaskedSpatialIntegralRAttrGG_
      end interface
      interface MaskedSpatialAverage ; module procedure &
	   MaskedSpatialAverageRAttrGG_
      end interface
      interface PairedSpatialIntegrals ; module procedure &
	    PairedSpatialIntegralRAttrGG_
      end interface
      interface PairedSpatialAverages ; module procedure &
	    PairedSpatialAverageRAttrGG_
      end interface
      interface PairedMaskedSpatialIntegrals ; module procedure &
	    PairedMaskedIntegralRAttrGG_
      end interface
      interface PairedMaskedSpatialAverages ; module procedure &
	    PairedMaskedAverageRAttrGG_
      end interface

! !REVISION HISTORY:
! 	25Oct01 - J.W. Larson <larson@mcs.anl.gov> - Initial version
!        9May02 - J.W. Larson <larson@mcs.anl.gov> - Massive Refactoring.
!    10-14Jun02 - J.W. Larson <larson@mcs.anl.gov> - Added Masked methods.
!    17-18Jun02 - J.W. Larson <larson@mcs.anl.gov> - Added Paired/Masked 
!                 methods.
!       18Jun02 - J.W. Larson <larson@mcs.anl.gov> - Renamed module from 
!                 m_GlobalIntegral to m_SpatialIntegral.
!       15Jan03 - E.T. Ong <eong@mcs.anl.gov> - Initialized real-only 
!                 AttrVects using nullfied integer lists. This circuitous 
!                 hack was required because the compaq compiler does not
!                 compile the function AttrVectExportListToChar. 
!
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname='MCT::m_SpatialIntegral'

 contains

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: SpatialIntegralRAttrGG_ - Compute spatial integral.
!
! !DESCRIPTION:
! This routine computes spatial integrals of the {\tt REAL} attributes
! of the {\tt REAL} attributes of the input {\tt AttrVect} argument 
! {\tt inAv}.  {\tt SpatialIntegralRAttrGG\_()} takes the input 
! {\tt AttrVect} argument {\tt inAv} and computes the spatial 
! integral using weights stored in the {\tt GeneralGrid} argument 
! {\tt GGrid} and identified by the {\tt CHARACTER} tag {\tt WeightTag}.  
! The integral of each {\tt REAL} attribute is returned in the output 
! {\tt AttrVect} argument {\tt outAv}. If {\tt SpatialIntegralRAttrGG\_()} 
! is invoked with the optional {\tt LOGICAL} input argument 
! {\tt SumWeights} set as {\tt .TRUE.}, then the weights are also summed 
! and stored in {\tt outAv} (and can be referenced with the attribute 
! tag defined by the argument{\tt WeightTag}.  If 
! {\tt SpatialIntegralRAttrGG\_()} is invoked with the optional {\tt INTEGER} 
! argument {\tt comm} (a Fortran MPI communicator handle), the summation
! operations for the integral are completed on the local process, then 
! reduced across the communicator, with all processes receiving the result.
!
! {\bf N.B.:  } The local lengths of the {\tt AttrVect} argument {\tt inAv} 
! and the {\tt GeneralGrid} {\tt GGrid} must be equal.  That is, there 
! must be a one-to-one correspondence between the field point values stored 
! in {\tt inAv} and the point weights stored in {\tt GGrid}.
!
! {\bf N.B.:  }  If {\tt SpatialIntegralRAttrGG\_()} is invoked with the 
! optional {\tt LOGICAL} input argument {\tt SumWeights} set as {\tt .TRUE.}, 
! then the value of {\tt WeightTag} must not conflict with any of the 
! {\tt REAL} attribute tags in {\tt inAv}.
!
! {\bf N.B.:  } The output {\tt AttrVect} argument {\tt outAv} is an 
! allocated data structure.  The user must deallocate it using the routine 
! {\tt AttrVect\_clean()} when it is no longer needed.  Failure to do so 
! will result in a memory leak.
!
! !INTERFACE:

 subroutine SpatialIntegralRAttrGG_(inAv, outAv, GGrid, WeightTag, &
                                   SumWeights, comm)
! ! USES:

      use m_stdio
      use m_die
      use m_mpif90

      use m_realkinds, only : FP

      use m_AttrVect, only : AttrVect
      use m_AttrVect, only : AttrVect_lsize => lsize

      use m_GeneralGrid, only : GeneralGrid
      use m_GeneralGrid, only : GeneralGrid_lsize => lsize
      use m_GeneralGrid, only : GeneralGrid_indexRA => indexRA
      use m_GeneralGrid, only : GeneralGrid_exportRAttr => exportRAttr

      use m_SpatialIntegralV, only: SpatialIntegralV

      implicit none

! !INPUT PARAMETERS:

      type(AttrVect),    intent(IN) :: inAv
      type(GeneralGrid), intent(IN) :: GGrid
      character(len=*),  intent(IN) :: WeightTag
      logical, optional, intent(IN) :: SumWeights
      integer, optional, intent(IN) :: comm

! !OUTPUT PARAMETERS:

      type(AttrVect),    intent(OUT) :: outAv

! !REVISION HISTORY:
! 	06Feb02 - J.W. Larson <larson@mcs.anl.gov> - initial version
! 	09May02 - J.W. Larson <larson@mcs.anl.gov> - Refactored and
!                 renamed SpatialIntegralRAttrGG_().
!       07Jun02 - J.W. Larson <larson@mcs.anl.gov> - Bug fix and further
!                 refactoring.
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::SpatialIntegralRAttrGG_'

  integer :: ierr, length
  logical :: mySumWeights
  real(FP), dimension(:), pointer :: gridWeights

       ! Argument Validity Checks

  if(AttrVect_lsize(inAv) /= GeneralGrid_lsize(GGrid)) then
     ierr = AttrVect_lsize(inAv) - GeneralGrid_lsize(GGrid)
     write(stderr,'(3a,i8,a,i8)') myname_, &
	  ':: inAv / GGrid length mismatch:  ', &
	  ' AttrVect_lsize(inAv) = ',AttrVect_lsize(inAv), &
	  ' GeneralGrid_lsize(GGrid) = ',GeneralGrid_lsize(GGrid)
     call die(myname_)
  endif

  if(present(SumWeights)) then
     mySumWeights = SumWeights
  else
     mySumWeights = .FALSE.
  endif

       ! ensure unambiguous pointer association status for gridWeights

  nullify(gridWeights) 

       ! Extract Grid Weights

  call GeneralGrid_exportRAttr(GGrid, WeightTag, gridWeights, length)

       ! 

  if(present(comm)) then ! do a distributed AllReduce-style integral:
     call SpatialIntegralV(inAv, outAv, gridWeights, mySumWeights, &
	                        WeightTag, comm)
  else
     call SpatialIntegralV(inAv, outAv, gridWeights, mySumWeights, &
                                WeightTag)
  endif

       ! Clean up temporary allocated space

  deallocate(gridWeights, stat=ierr)
  if(ierr /= 0) then
     write(stderr,'(2a,i8)') myname_, &
	              ':: deallocate(gridWeights...failed.  ierr=', ierr
     call die(myname_)
  endif

 end subroutine SpatialIntegralRAttrGG_

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: SpatialAverageRAttrGG_ - Compute spatial average.
!
! !DESCRIPTION:
! This routine computes spatial averages of the {\tt REAL} attributes
! of the input {\tt AttrVect} argument {\tt inAv}.  
! {\tt SpatialAverageRAttrGG\_()} takes the input {\tt AttrVect} argument 
! {\tt inAv} and computes the spatial average using weights 
! stored in the {\tt GeneralGrid} argument {\tt GGrid} and identified by
! the {\tt CHARACTER} tag {\tt WeightTag}.  The average of each {\tt REAL}
! attribute is returned in the output {\tt AttrVect} argument {\tt outAv}.
! If {\tt SpatialAverageRAttrGG\_()} is invoked with the optional {\tt INTEGER} 
! argument {\tt comm} (a Fortran MPI communicator handle), the summation
! operations for the average are completed on the local process, then 
! reduced across the communicator, with all processes receiving the result.
!
! {\bf N.B.:  } The local lengths of the {\tt AttrVect} argument {\tt inAv} 
! and the {\tt GeneralGrid} {\tt GGrid} must be equal.  That is, there 
! must be a one-to-one correspondence between the field point values stored 
! in {\tt inAv} and the point weights stored in {\tt GGrid}.
!
! {\bf N.B.:  } The output {\tt AttrVect} argument {\tt outAv} is an 
! allocated data structure.  The user must deallocate it using the routine 
! {\tt AttrVect\_clean()} when it is no longer needed.  Failure to do so 
! will result in a memory leak.
!
! !INTERFACE:

 subroutine SpatialAverageRAttrGG_(inAv, outAv, GGrid, WeightTag, comm)

! ! USES:

      use m_realkinds, only : FP
   
      use m_stdio
      use m_die
      use m_mpif90

      use m_AttrVect, only : AttrVect
      use m_AttrVect, only : AttrVect_init => init
      use m_AttrVect, only : AttrVect_zero => zero 
      use m_AttrVect, only : AttrVect_clean => clean
      use m_AttrVect, only : AttrVect_nRAttr => nRAttr
      use m_AttrVect, only : AttrVect_indexRA => indexRA

      use m_GeneralGrid, only : GeneralGrid

      use m_List, only : List
      use m_List, only : List_nullify => nullify

      implicit none

! !INPUT PARAMETERS:

      type(AttrVect),    intent(IN) :: inAv
      type(GeneralGrid), intent(IN) :: GGrid
      character(len=*),  intent(IN) :: WeightTag
      integer, optional, intent(IN) :: comm

! !OUTPUT PARAMETERS:

      type(AttrVect),    intent(OUT) :: outAv

! !REVISION HISTORY:
! 	08Feb02 - J.W. Larson <larson@mcs.anl.gov> - initial version
! 	08May02 - J.W. Larson <larson@mcs.anl.gov> - minor modifications:
!                 1) renamed the routine to GlobalAverageRAttrGG_
!                 2) changed calls to reflect new routine name
!                    GlobalIntegralRAttrGG_().
!       18Jun02 - J.W. Larson <larson@mcs.anl.gov> - Renamed routine to
!                 SpatialAverageRAttrGG_().
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::SpatialAverageRAtttrGG_'

  type(AttrVect) :: integratedAv
  type(List) :: nullIList
  integer :: i, ierr, iweight

       ! Compute the spatial integral:

  if(present(comm)) then
     call SpatialIntegralRAttrGG_(inAv, integratedAv, GGrid, WeightTag, &
	                         .TRUE., comm)
  else
     call SpatialIntegralRAttrGG_(inAv, integratedAv, GGrid, WeightTag, &
	                         .TRUE.)
  endif

       ! Check value of summed weights (to avoid division by zero):

  iweight = AttrVect_indexRA(integratedAv, WeightTag)
  if(integratedAv%rAttr(iweight, 1) == 0._FP) then
     write(stderr,'(2a)') myname_, &
	  '::ERROR--Global sum of grid weights is zero.'
     call die(myname_)
  endif

       ! Initialize output AttrVect outAv:

  call List_nullify(nullIList)
  call AttrVect_init(outAv, iList=nullIList, rList=inAv%rList, lsize=1)
  call AttrVect_zero(outAv)

       ! Divide by global weight sum to compute spatial averages from 
       ! spatial integrals.

  do i=1,AttrVect_nRAttr(outAv)
     outAv%rAttr(i,1) = integratedAv%rAttr(i,1) &
	                               / integratedAv%rAttr(iweight,1) 
  end do

       ! Clean up temporary AttrVect:

  call AttrVect_clean(integratedAv)

 end subroutine SpatialAverageRAttrGG_

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: MaskedSpatialIntegralRAttrGG_ - Masked spatial integral.
!
! !DESCRIPTION: 
! This routine computes masked spatial integrals of the {\tt REAL} 
! attributes of the input {\tt AttrVect} argument {\tt inAv}, returning 
! the masked integrals in the output {\tt AttrVect} {\tt outAv}.  All of 
! the masking data are assumed stored in the input {\tt GeneralGrid} 
! argument {\tt GGrid}.  If integer masks are to be used, their integer 
! attribute names in {\tt GGrid} are named as a colon-delimited list 
! in the optional {\tt CHARACTER} input argument {\tt iMaskTags}.  Real 
! masks (if desired) are referenced by their real attribute names in 
! {\tt GGrid} are named as a colon-delimited list in the optional 
! {\tt CHARACTER} input argument {\tt rMaskTags}.  The user specifies 
! a choice of mask combination method with the input {\tt LOGICAL} argument
! {\tt UseFastMethod}.  If ${\tt UseFastMethod} = {\tt .FALSE.}$ this 
! routine checks each mask entry to ensure that the integer masks contain
! only ones and zeroes, and that entries in the real masks are all in
! the closed interval $[0,1]$.  If ${\tt UseFastMethod} = {\tt .TRUE.}$,
! this routine performs direct products of the masks, assuming that the
! user has validated them in advance.  The optional {\tt LOGICAL} input 
! argument {\tt SumWeights} determines whether the masked sum of the spatial
! weights is computed and returned in {\tt outAv} with the real attribute 
! name supplied in the optional {\tt CHARACTER} input argument 
! {\tt WeightSumTag}.  This integral can either be a local (i.e. a global 
! memory space operation), or a global distributed integral.  The latter 
! is the case if the optional input {\tt INTEGER} argument {\tt comm} is
! supplied (which corresponds to a Fortran MPI communicatior handle).
!
! {\bf N.B.:  } The local lengths of the {\tt AttrVect} argument {\tt inAv} 
! and the input {\tt GeneralGrid} {\tt GGrid} must be equal.  That is, there 
! must be a one-to-one correspondence between the field point values stored 
! in {\tt inAv} and the point weights stored in {\tt GGrid}.
!
! {\bf N.B.:  }  If {\tt SpatialIntegralRAttrV\_()} is invoked with the 
! optional {\tt LOGICAL} input argument {\tt SumWeights} set as {\tt .TRUE.}. 
! In this case, the none of {\tt REAL} attribute tags in {\tt inAv} may be 
! named the same as the string contained in {\tt WeightSumTag}, which is an 
! attribute name reserved for the sum of the weights in the output {\tt AttrVect} 
! {\tt outAv}.
!
! {\bf N.B.:  } The output {\tt AttrVect} argument {\tt outAv} is an 
! allocated data structure.  The user must deallocate it using the routine 
! {\tt AttrVect\_clean()} when it is no longer needed.  Failure to do so 
! will result in a memory leak.
!
! !INTERFACE:

 subroutine MaskedSpatialIntegralRAttrGG_(inAv, outAv, GGrid, SpatialWeightTag, &
                                         iMaskTags, rMaskTags, UseFastMethod,  &
                                         SumWeights, WeightSumTag, comm)

! ! USES:

      use m_stdio
      use m_die
      use m_mpif90

      use m_realkinds, only : FP

      use m_String, only : String
      use m_String, only : String_toChar => toChar
      use m_String, only : String_clean => clean

      use m_List, only : List
      use m_List, only : List_init => init
      use m_List, only : List_clean => clean
      use m_List, only : List_nitem => nitem
      use m_List, only : List_get => get

      use m_AttrVect, only : AttrVect
      use m_AttrVect, only : AttrVect_lsize => lsize

      use m_GeneralGrid, only : GeneralGrid
      use m_GeneralGrid, only : GeneralGrid_lsize => lsize
      use m_GeneralGrid, only : GeneralGrid_indexRA => indexRA
      use m_GeneralGrid, only : GeneralGrid_exportIAttr => exportIAttr
      use m_GeneralGrid, only : GeneralGrid_exportRAttr => exportRAttr

      use m_AttrVectReduce, only : AttrVect_GlobalWeightedSumRAttr => &
	                                         GlobalWeightedSumRAttr
      use m_AttrVectReduce, only : AttrVect_LocalWeightedSumRAttr => &
	                                         LocalWeightedSumRAttr

      use m_SpatialIntegralV, only : MaskedSpatialIntegralV

      implicit none

! !INPUT PARAMETERS:

      type(AttrVect),                  intent(IN) :: inAv
      type(GeneralGrid),               intent(IN) :: GGrid
      character(len=*),                intent(IN) :: SpatialWeightTag
      character(len=*),      optional, intent(IN) :: iMaskTags
      character(len=*),      optional, intent(IN) :: rMaskTags
      logical,                         intent(IN) :: UseFastMethod
      logical,               optional, intent(IN) :: SumWeights
      character(len=*),      optional, intent(IN) :: WeightSumTag
      integer,               optional, intent(IN) :: comm

! !OUTPUT PARAMETERS:

      type(AttrVect),                  intent(OUT) :: outAv

! !REVISION HISTORY:
! 	11Jun02 - J.W. Larson <larson@mcs.anl.gov> - initial version
!
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::MaskedSpatialIntegralRAttrGG_'

  integer :: i, ierr, j, length
  logical :: mySumWeights

  type(List) :: iMaskList, rMaskList
  type(String) :: DummStr

  integer, dimension(:), pointer :: iMask, iMaskTemp
  real(FP), dimension(:), pointer :: rMask, rMaskTemp
  integer :: TempMaskLength

  real(FP), dimension(:), pointer :: SpatialWeights

  integer :: niM, nrM ! Number of iMasks and rMasks, respectively

       ! Argument Validity Checks

  if(AttrVect_lsize(inAv) /= GeneralGrid_lsize(GGrid)) then
     ierr = AttrVect_lsize(inAv) - GeneralGrid_lsize(GGrid)
     write(stderr,'(3a,i8,a,i8)') myname_, &
	  ':: inAv / GGrid length mismatch:  ', &
	  ' AttrVect_lsize(inAv) = ',AttrVect_lsize(inAv), &
	  ' GeneralGrid_lsize(GGrid) = ',GeneralGrid_lsize(GGrid)
     call die(myname_)
  endif

  if(present(SumWeights)) then
     mySumWeights = SumWeights
     if(.not. present(WeightSumTag)) then
	write(stderr,'(3a)') myname_,':: FATAL--If the input argument SumWeights=.TRUE.,', &
                        ' then the argument WeightSumTag must be provided.'
	call die(myname_)
     endif
  else
     mySumWeights = .FALSE.
  endif

  if(present(iMaskTags)) then
     call List_init(iMaskList, iMaskTags)
     if(List_nitem(iMaskList) == 0) then
	write(stderr,'(3a)') myname_,':: ERROR--an INTEGER mask list with', &
	               'no valid items was provided.'
        call die(myname_)
     endif
  endif

  if(present(rMaskTags)) then
     call List_init(rMaskList, rMaskTags)
     if(List_nitem(iMaskList) == 0) then
	write(stderr,'(3a)') myname_,':: ERROR--an REAL mask list with', &
	               'no valid items was provided.'
        call die(myname_)
     endif
  endif

       ! Determine the on-processor vector length for use throughout
       ! this routine:

  length = AttrVect_lsize(inAv)

       !==========================================================
       ! Extract Spatial Weights from GGrid using SpatialWeightTag
       !==========================================================

  nullify(SpatialWeights)
  call GeneralGrid_exportRAttr(GGrid, SpatialWeightTag, SpatialWeights, &
                               TempMaskLength)
  if(TempMaskLength /= length) then
     write(stderr,'(3a,i8,a,i8)') myname_,&
	  ':: error on return from GeneralGrid_exportRAttr().' , &   
	  'Returned with SpatialWeights(:) length = ',TempMaskLength, &
	  ',which conflicts with AttrVect_lsize(inAv) = ',length
     call die(myname_)
  endif

       !==========================================================
       ! If the argument iMaskTags is present, create the combined
       ! iMask array:
       !==========================================================

  if(present(iMaskTags)) then ! assemble iMask(:) from all the integer
                              ! mask attributes stored in GGrid(:)

     allocate(iMask(length), iMaskTemp(length), stat=ierr)
     if(ierr /= 0) then
	write(stderr,'(3a,i8)') myname_,':: allocate(iMask(...) failed,', &
	     ' ierr=',ierr
	call die(myname_)
     endif

     niM = List_nitem(iMaskList)

     do i=1,niM

       ! Retrieve current iMask tag, and get this attribute from GGrid:
	call List_get(DummStr, i, iMaskList)
	call GeneralGrid_exportIAttr(GGrid, String_toChar(DummStr), &
	                             iMaskTemp, TempMaskLength)
	call String_clean(DummStr)
	if(TempMaskLength /= length) then
	   write(stderr,'(3a,i8,a,i8)') myname_,&
		 ':: error on return from GeneralGrid_exportIAttr().' , &   
		 'Returned with TempMaskLength = ',TempMaskLength, &
		 ',which conflicts with AttrVect_lsize(inAv) = ',length
	   call die(myname_)
	endif

	if(i == 1) then ! first pass--examine iMaskTemp(:) only

	   if(UseFastMethod) then ! straight copy of iMaskTemp(:)
	      do j=1,length
		 iMask(j) = iMaskTemp(j)
	      end do
	   else ! go through the entries of iMaskTemp(:) one-by-one
	      do j=1,length
		 select case(iMaskTemp(j))
		 case(0)
		    iMask(j) = 0
		 case(1)
		    iMask(j) = 1
		 case default
		    write(stderr,'(3a,i8,a,i8)') myname_, &
			 ':: FATAL--illegal INTEGER mask entry.  Integer mask ', &
			 'entries must be 0 or 1.  iMask(',j,') = ', iMask(j)
		    call die(myname_)
		 end select ! select case(iMaskTemp(j))...
	      end do ! do j=1,length
	   endif ! if(UseFastMethod)...

	else ! That is, i /= 1 ...

	   if(UseFastMethod) then ! straight product of iMask(:) 
                                  ! and iMaskTemp(:)
	      do j=1,length
		 iMask(j) = iMask(j) * iMaskTemp(j)
	      end do
	   else ! go through the entries of iMaskTemp(:) one-by-one
	      do j=1,length
		 select case(iMaskTemp(j))
		 case(0) ! zero out iMask(j)
		    iMask(j) = 0
		 case(1) ! do nothing
		 case default
		    write(stderr,'(3a,i8,a,i8)') myname_, &
			 ':: FATAL--illegal INTEGER mask entry.  Integer mask ', &
			 'entries must be 0 or 1.  iMask(',j,') = ', iMask(j)
		    call die(myname_)
		 end select ! select case(iMaskTemp(j))...
	      end do ! do j=1,length
	   endif ! if(UseFastMethod)...

	endif ! if(i == 1)...

     end do ! do i=1,niM...iMask retrievals

  endif ! if(present(iMaskTags))...

       !==========================================================
       ! If the argument rMaskTags is present, create the combined
       ! REAL mask rMask array:
       !==========================================================

  if(present(rMaskTags)) then ! assemble rMask(:) from all the integer
                              ! mask attributes stored in GGrid(:)

     allocate(rMask(length), rMaskTemp(length), stat=ierr)
     if(ierr /= 0) then
	write(stderr,'(3a,i8)') myname_,':: allocate(rMask(...) failed,', &
	     ' ierr=',ierr
	call die(myname_)
     endif

     nrM = List_nitem(rMaskList)

     do i=1,nrM

       ! Retrieve current rMask tag, and get this attribute from GGrid:
	call List_get(DummStr, i, rMaskList)
	call GeneralGrid_exportRAttr(GGrid, String_toChar(DummStr), &
	                             rMaskTemp, TempMaskLength)
	call String_clean(DummStr)
	if(TempMaskLength /= length) then
	   write(stderr,'(3a,i8,a,i8)') myname_,&
		 ':: error on return from GeneralGrid_exportRAttr().' , &   
		 'Returned with TempMaskLength = ',TempMaskLength, &
		 ',which conflicts with AttrVect_lsize(inAv) = ',length
	   call die(myname_)
	endif

	if(i == 1) then ! first pass--examine rMaskTemp(:) only

	   if(UseFastMethod) then ! straight copy of rMaskTemp(:)
	      do j=1,length
		 rMask(j) = rMaskTemp(j)
	      end do
	   else ! go through the entries of rMaskTemp(:) one-by-one
                ! to ensure they are in the range [0.,1.]
	      do j=1,length
		 if((rMaskTemp(j) >= 0.) .or. (rMaskTemp(j) <=1.)) then
		    rMask(j) = rMaskTemp(j)
		 else
		    write(stderr,'(3a,i8,a,i8)') myname_, &
			 ':: FATAL--illegal REAL mask entry.  Real mask ', &
			 'entries must be in [0.,1.]  rMask(',j,') = ', rMask(j)
		    call die(myname_)
		 endif ! if((rMaskTemp(j) >= 0.) .or. (rMaskTemp(j) <=1.))...
	      end do ! do j=1,length
	   endif ! if(UseFastMethod)...

	else ! That is, i /= 1 ...

	   if(UseFastMethod) then ! straight product of rMask(:) 
                                  ! and rMaskTemp(:)
	      do j=1,length
		 rMask(j) = rMask(j) * rMaskTemp(j)
	      end do
	   else ! go through the entries of rMaskTemp(:) one-by-one
                ! to ensure they are in the range [0.,1.]
	      do j=1,length
		 if((rMaskTemp(j) >= 0.) .or. (rMaskTemp(j) <=1.)) then
		    rMask(j) = rMask(j) * rMaskTemp(j)
		 else
		    write(stderr,'(3a,i8,a,i8)') myname_, &
			 ':: FATAL--illegal REAL mask entry.  Real mask ', &
			 'entries must be in [0.,1.]  rMask(',j,') = ', rMask(j)
		    call die(myname_)
		 endif ! if((rMaskTemp(j) >= 0.) .or. (rMaskTemp(j) <=1.))...
	      end do ! do j=1,length
	   endif ! if(UseFastMethod)...

	endif ! if(i == 1)...

     end do ! do i=1,niM...rMask retrievals

  endif ! if(present(rMaskTags))...

       !==========================================================
       ! Now that we have produced single INTEGER and REAL masks,
       ! compute the masked weighted sum.
       !==========================================================

  if(present(rMaskTags)) then ! We have a REAL Mask

     if(present(iMaskTags)) then ! and an INTEGER Mask

	if(present(comm)) then ! compute distributed AllReduce-style sum:
	   
	   if(mySumWeights) then ! return the global masked sum of the 
                                 ! weights in outAV
	      call MaskedSpatialIntegralV(inAv, outAv, SpatialWeights, &
		                               iMask, rMask, UseFastMethod, &
					       SumWeights, WeightSumTag, comm)
	   else ! Do not return the masked sum of the weights
	      call MaskedSpatialIntegralV(inAv, outAv, SpatialWeights, &
		                               iMask, rMask, UseFastMethod, &
					       comm=comm)
	   endif ! if(mySumWeights)...

	else ! compute local sum:

	   if(mySumWeights) then ! return the global masked sum of the 
                                 ! weights in outAV
	      call MaskedSpatialIntegralV(inAv, outAv, SpatialWeights, &
		                               iMask, rMask, UseFastMethod, &
					       SumWeights, WeightSumTag)
	   else ! Do not return the masked sum of the weights
	      call MaskedSpatialIntegralV(inAv, outAv, SpatialWeights, &
		                               iMask, rMask, UseFastMethod)
	   endif ! if(mySumWeights)...

	endif ! if(present(comm))...

     else ! REAL Mask Only Case...

	if(present(comm)) then ! compute distributed AllReduce-style sum:
	   
	   if(mySumWeights) then ! return the global masked sum of the 
                                 ! weights in outAV
	      call MaskedSpatialIntegralV(inAv, outAv, SpatialWeights, &
		                               rMask=rMask, &
					       UseFastMethod=UseFastMethod, &
					       SumWeights=SumWeights, &
					       WeightSumTag=WeightSumTag, &
					       comm=comm)
	   else ! Do not return the masked sum of the weights
	      call MaskedSpatialIntegralV(inAv, outAv, SpatialWeights, &
		                               rMask=rMask, &
					       UseFastMethod=UseFastMethod, &
					       comm=comm)
	   endif ! if(mySumWeights)...

	else ! compute local sum:

	   if(mySumWeights) then ! return the global masked sum of the 
                                 ! weights in outAV
	      call MaskedSpatialIntegralV(inAv, outAv, SpatialWeights, &
		                               rMask=rMask, &
					       UseFastMethod=UseFastMethod, &
					       SumWeights=SumWeights, &
					       WeightSumTag=WeightSumTag)
	   else ! Do not return the masked sum of the weights
	      call MaskedSpatialIntegralV(inAv, outAv, SpatialWeights, &
		                               rMask=rMask, &
					       UseFastMethod=UseFastMethod)
	   endif ! if(mySumWeights)...

	endif ! if(present(comm))...

     endif
  else ! no REAL Mask...

     if(present(iMaskTags)) then ! INTEGER Mask Only Case...

	if(present(comm)) then ! compute distributed AllReduce-style sum:
	   
	   if(mySumWeights) then ! return the global masked sum of the 
                                 ! weights in outAV
	      call MaskedSpatialIntegralV(inAv, outAv, SpatialWeights, &
		                               iMask=iMask, &
					       UseFastMethod=UseFastMethod, &
					       SumWeights=SumWeights, &
					       WeightSumTag=WeightSumTag, &
					       comm=comm)
	   else ! Do not return the masked sum of the weights
	      call MaskedSpatialIntegralV(inAv, outAv, SpatialWeights, &
		                               iMask=iMask, &
					       UseFastMethod=UseFastMethod, &
					       comm=comm)
	   endif ! if(mySumWeights)...

	else ! compute local sum:

	   if(mySumWeights) then ! return the global masked sum of the 
                                 ! weights in outAV
	      call MaskedSpatialIntegralV(inAv, outAv, SpatialWeights, &
		                               iMask=iMask, &
					       UseFastMethod=UseFastMethod, &
					       SumWeights=SumWeights, &
					       WeightSumTag=WeightSumTag)
	   else ! Do not return the masked sum of the weights
	      call MaskedSpatialIntegralV(inAv, outAv, SpatialWeights, &
		                               iMask=iMask, &
					       UseFastMethod=UseFastMethod)
	   endif ! if(mySumWeights)...

	endif ! if(present(comm))...

     else ! no INTEGER Mask / no REAL Mask Case...

	if(present(comm)) then ! compute distributed AllReduce-style sum:
	   
	   if(mySumWeights) then ! return the global masked sum of the 
                                 ! weights in outAV
	      call MaskedSpatialIntegralV(inAv, outAv, SpatialWeights, &
					       UseFastMethod=UseFastMethod, &
					       SumWeights=SumWeights, &
					       WeightSumTag=WeightSumTag, &
					       comm=comm)
	   else ! Do not return the masked sum of the weights
	      call MaskedSpatialIntegralV(inAv, outAv, SpatialWeights, &
					       UseFastMethod=UseFastMethod, &
					       comm=comm)
	   endif ! if(mySumWeights)...

	else ! compute local sum:

	   if(mySumWeights) then ! return the global masked sum of the 
                                 ! weights in outAV
	      call MaskedSpatialIntegralV(inAv, outAv, SpatialWeights, &
					       UseFastMethod=UseFastMethod, &
					       SumWeights=SumWeights, &
					       WeightSumTag=WeightSumTag)
	   else ! Do not return the masked sum of the weights
	      call MaskedSpatialIntegralV(inAv, outAv, SpatialWeights, &
					       UseFastMethod=UseFastMethod)
	   endif ! if(mySumWeights)...

	endif ! if(present(comm))...

     endif ! if(present(iMaskTags)...

  endif ! if(present(rMaskTags)...

       !==========================================================
       ! The masked spatial integral is now completed.
       ! Clean up the the various allocated mask structures.
       !==========================================================

  if(present(iMaskTags)) then ! clean up iMask and friends...
     call List_clean(iMaskList)
     deallocate(iMask, iMaskTemp, stat=ierr)
     if(ierr /= 0) then
	write(stderr,'(3a,i8)') myname_,':: deallocate(iMask(...) failed,', &
	     ' ierr=',ierr
	call die(myname_)
     endif
  endif

  if(present(rMaskTags)) then ! clean up rMask and co...
     call List_clean(rMaskList)
     deallocate(rMask, rMaskTemp, stat=ierr)
     if(ierr /= 0) then
	write(stderr,'(3a,i8)') myname_,':: deallocate(rMask(...) failed,', &
	     ' ierr=',ierr
	call die(myname_)
     endif
  endif

       ! Clean up SpatialWeights(:)

  deallocate(SpatialWeights, stat=ierr)
  if(ierr /= 0) then
     write(stderr,'(3a,i8)') myname_,':: deallocate(SpatialWeights(...) failed,', &
	                    ' ierr=',ierr
     call die(myname_)
  endif

 end subroutine MaskedSpatialIntegralRAttrGG_

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: MaskedSpatialAverageRAttrGG_ - Masked spatial average.
!
! !DESCRIPTION: 
! This routine computes masked spatial averages of the {\tt REAL} 
! attributes of the input {\tt AttrVect} argument {\tt inAv}, returning 
! the masked averages in the output {\tt AttrVect} {\tt outAv}.  All of 
! the masking data are assumed stored in the input {\tt GeneralGrid} 
! argument {\tt GGrid}.  If integer masks are to be used, their integer 
! attribute names in {\tt GGrid} are named as a colon-delimited list 
! in the optional {\tt CHARACTER} input argument {\tt iMaskTags}.  Real 
! masks (if desired) are referenced by their real attribute names in 
! {\tt GGrid} are named as a colon-delimited list in the optional 
! {\tt CHARACTER} input argument {\tt rMaskTags}.  The user specifies 
! a choice of mask combination method with the input {\tt LOGICAL} argument
! {\tt UseFastMethod}.  If ${\tt UseFastMethod} = {\tt .FALSE.}$ this 
! routine checks each mask entry to ensure that the integer masks contain
! only ones and zeroes, and that entries in the real masks are all in
! the closed interval $[0,1]$.  If ${\tt UseFastMethod} = {\tt .TRUE.}$,
! this routine performs direct products of the masks, assuming that the
! user has validated them in advance.  This averaging can either be a 
! local (equivalent to a global memory space operation), or a global 
! distributed integral.  The latter is the case if the optional input 
! {\tt INTEGER} argument {\tt comm} is supplied (which corresponds to a 
! Fortran MPI communicatior handle).
!
! {\bf N.B.:  } The local lengths of the {\tt AttrVect} argument {\tt inAv} 
! and the input {\tt GeneralGrid} {\tt GGrid} must be equal.  That is, 
! there must be a one-to-one correspondence between the field point values 
! stored in {\tt inAv} and the point weights stored in {\tt GGrid}.
!
! {\bf N.B.:  } The output {\tt AttrVect} argument {\tt outAv} is an 
! allocated data structure.  The user must deallocate it using the routine 
! {\tt AttrVect\_clean()} when it is no longer needed.  Failure to do so 
! will result in a memory leak.
!
! !INTERFACE:

 subroutine MaskedSpatialAverageRAttrGG_(inAv, outAv, GGrid, SpatialWeightTag, &
                                        iMaskTags, rMaskTags, UseFastMethod,  &
  				        comm)

! ! USES:

      use m_realkinds, only : FP

      use m_stdio
      use m_die
      use m_mpif90

      use m_AttrVect, only : AttrVect
      use m_AttrVect, only : AttrVect_init => init
      use m_AttrVect, only : AttrVect_zero => zero 
      use m_AttrVect, only : AttrVect_clean => clean
      use m_AttrVect, only : AttrVect_lsize => lsize
      use m_AttrVect, only : AttrVect_indexRA => indexRA
      use m_AttrVect, only : AttrVect_nRAttr => nRAttr

      use m_GeneralGrid, only : GeneralGrid
      use m_GeneralGrid, only : GeneralGrid_lsize => lsize
      use m_GeneralGrid, only : GeneralGrid_indexRA => indexRA

      use m_List, only : List
      use m_List, only : List_nullify => nullify

      implicit none

! !INPUT PARAMETERS:

      type(AttrVect),                  intent(IN) :: inAv
      type(GeneralGrid),               intent(IN) :: GGrid
      character(len=*),                intent(IN) :: SpatialWeightTag
      character(len=*),      optional, intent(IN) :: iMaskTags
      character(len=*),      optional, intent(IN) :: rMaskTags
      logical,                         intent(IN) :: UseFastMethod
      integer,               optional, intent(IN) :: comm

! !OUTPUT PARAMETERS:

      type(AttrVect),                  intent(OUT) :: outAv

! !REVISION HISTORY:
! 	12Jun02 - J.W. Larson <larson@mcs.anl.gov> - initial version
!
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::MaskedSpatialAverageRAttrGG_'

  type(AttrVect) :: integratedAv
  type(List) :: nullIList
  character*9, parameter :: WeightSumTag = 'WeightSum'

  integer :: i, iweight

       !================================================================
       ! Do the integration using MaskedSpatialIntegralRAttrGG_(), which
       ! returns the intermediate integrals (including the masked weight
       ! sum) in the AttrVect integratedAv.
       !================================================================

  if(present(iMaskTags)) then

     if(present(rMaskTags)) then ! have both iMasks and rMasks

	if(present(comm)) then ! a distributed parallel sum
	   call MaskedSpatialIntegralRAttrGG_(inAv, integratedAv, GGrid, &
		                             SpatialWeightTag, iMaskTags, &
					     rMaskTags, UseFastMethod,  &
					     .TRUE., WeightSumTag, comm)
	else ! a purely local sum
	   call MaskedSpatialIntegralRAttrGG_(inAv, integratedAv, GGrid, &
		                             SpatialWeightTag, iMaskTags, &
					     rMaskTags, UseFastMethod,  &
					     .TRUE., WeightSumTag)
	endif ! if(present(comm))...

     else ! Only iMasks are in use

	if(present(comm)) then ! a distributed parallel sum
	   call MaskedSpatialIntegralRAttrGG_(inAv, integratedAv, GGrid, &
		                             SpatialWeightTag, iMaskTags, &
					     UseFastMethod=UseFastMethod,  &
					     SumWeights=.TRUE., &
					     WeightSumTag=WeightSumTag, &
					     comm=comm)

	else ! a purely local sum
	   call MaskedSpatialIntegralRAttrGG_(inAv, integratedAv, GGrid, &
		                             SpatialWeightTag, iMaskTags, &
					     UseFastMethod=UseFastMethod,  &
					     SumWeights=.TRUE., &
					     WeightSumTag=WeightSumTag)
	endif ! if(present(comm))...

     endif ! if(present(rMaskTags)...

  else ! no iMasks

     if(present(rMaskTags)) then ! Only rMasks are in use

	if(present(comm)) then ! a distributed parallel sum
	   call MaskedSpatialIntegralRAttrGG_(inAv, integratedAv, GGrid, &
		                             SpatialWeightTag, &
					     rMaskTags=rMaskTags, &
					     UseFastMethod=UseFastMethod,  &
					     SumWeights=.TRUE., &
					     WeightSumTag=WeightSumTag, &
					     comm=comm)
	else ! a purely local sum
	   call MaskedSpatialIntegralRAttrGG_(inAv, integratedAv, GGrid, &
		                             SpatialWeightTag, &
					     rMaskTags=rMaskTags, &
					     UseFastMethod=UseFastMethod,  &
					     SumWeights=.TRUE., &
					     WeightSumTag=WeightSumTag)
	endif

     else ! Neither iMasks nor rMasks are in use

	if(present(comm)) then ! a distributed parallel sum
	   call MaskedSpatialIntegralRAttrGG_(inAv, integratedAv, GGrid, &
		                             SpatialWeightTag, &
					     UseFastMethod=UseFastMethod,  &
					     SumWeights=.TRUE., &
					     WeightSumTag=WeightSumTag, &
					     comm=comm)
	else ! a purely local sum
	   call MaskedSpatialIntegralRAttrGG_(inAv, integratedAv, GGrid, &
		                             SpatialWeightTag, &
					     UseFastMethod=UseFastMethod,  &
					     SumWeights=.TRUE., &
					     WeightSumTag=WeightSumTag)
	endif ! if(present(comm))...

     endif ! if(present(rMaskTags))...

  endif ! if(present(iMaskTags))...

       !================================================================
       ! The masked integrals and masked weight sum now reside in 
       ! in the AttrVect integratedAv.  We now wish to compute the
       ! averages by dividing the integtrals by the masked weight sum.
       !================================================================

       ! Check value of summed weights (to avoid division by zero):

  iweight = AttrVect_indexRA(integratedAv, WeightSumTag)
  if(integratedAv%rAttr(iweight, 1) == 0._FP) then
     write(stderr,'(2a)') myname_, &
	  '::ERROR--Global sum of grid weights is zero.'
     call die(myname_)
  endif

       ! Initialize output AttrVect outAv:
  call List_nullify(nullIList)
  call AttrVect_init(outAv, iList=nullIList, rList=inAv%rList, lsize=1)
  call AttrVect_zero(outAv)

       ! Divide by global weight sum to compute spatial averages from 
       ! spatial integrals.

  do i=1,AttrVect_nRAttr(outAv)
     outAv%rAttr(i,1) = integratedAv%rAttr(i,1) &
	                               / integratedAv%rAttr(iweight,1) 
  end do

       ! Clean up temporary AttrVect:

  call AttrVect_clean(integratedAv)

 end subroutine MaskedSpatialAverageRAttrGG_


!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: PairedSpatialIntegralRAttrGG_ - Do two spatial integrals at once.
!
! !DESCRIPTION:
! This routine computes spatial integrals of the {\tt REAL} attributes
! of the {\tt REAL} attributes of the input {\tt AttrVect} arguments 
! {\tt inAv1} and {\tt inAv2}, returning the integrals in the output 
! {\tt AttrVect} arguments {\tt outAv1} and {\tt outAv2}, respectively .  
! The integrals of {\tt inAv1} and {\tt inAv2} are computed using 
! spatial weights stored in the input {\tt GeneralGrid} arguments  
! {\tt GGrid1} and {\tt GGrid2}, respectively.  The spatial weights in 
! in {\tt GGrid1} and {\tt GGrid2} are identified by the input {\tt CHARACTER} 
! arguments {\tt WeightTag1} and {\tt WeightTag2}, respectively.  
! If {\tt SpatialIntegralRAttrGG\_()} is invoked with the optional 
! {\tt LOGICAL} input argument 
! {\tt SumWeights} set as {\tt .TRUE.}, then the weights are also summed 
! and stored in {\tt outAv1} and {\tt outAv2}, and can be referenced with 
! the attribute tags defined by the arguments {\tt WeightTag1} and 
! {\tt WeightTag2}, respectively.  This paired integral is implicitly a 
! distributed operation (the whole motivation for pairing the integrals is 
! to reduce communication latency costs), and the Fortran MPI communicator
! handle is defined by the input {\tt INTEGER} argument {\tt comm}.  The 
! summation is an AllReduce operation, with all processes receiving the 
! global sum.
!
! {\bf N.B.:  } The local lengths of the {\tt AttrVect} argument {\tt inAv1} 
! and the {\tt GeneralGrid} {\tt GGrid1} must be equal.  That is, there 
! must be a one-to-one correspondence between the field point values stored 
! in {\tt inAv1} and the point weights stored in {\tt GGrid1}.  The same
! relationship must apply between {\tt inAv2} and {\tt GGrid2}.
!
! {\bf N.B.:  }  If {\tt SpatialIntegralRAttrGG\_()} is invoked with the 
! optional {\tt LOGICAL} input argument {\tt SumWeights} set as {\tt .TRUE.}, 
! then the value of {\tt WeightTag1} must not conflict with any of the 
! {\tt REAL} attribute tags in {\tt inAv1} and the value of {\tt WeightTag2} 
! must not conflict with any of the {\tt REAL} attribute tags in {\tt inAv2}.
!
! {\bf N.B.:  } The output {\tt AttrVect} arguments {\tt outAv1} and 
! {\tt outAv2} are allocated data structures.  The user must deallocate them
!  using the routine {\tt AttrVect\_clean()} when they are no longer needed.  
! Failure to do so will result in a memory leak.
!
! !INTERFACE:

 subroutine PairedSpatialIntegralRAttrGG_(inAv1, outAv1, GGrid1, WeightTag1, &
                                         inAv2, outAv2, GGrid2, WeightTag2, &
					 SumWeights, comm)
! ! USES:

      use m_stdio
      use m_die
      use m_mpif90

      use m_realkinds, only : FP

      use m_AttrVect, only : AttrVect
      use m_AttrVect, only : AttrVect_lsize => lsize
      use m_AttrVect, only : AttrVect_nRAttr => nRAttr

      use m_GeneralGrid, only : GeneralGrid
      use m_GeneralGrid, only : GeneralGrid_lsize => lsize
      use m_GeneralGrid, only : GeneralGrid_indexRA => indexRA
      use m_GeneralGrid, only : GeneralGrid_exportRAttr => exportRAttr

      use m_AttrVectReduce, only : AttrVect_LocalWeightedSumRAttr => &
	                                         LocalWeightedSumRAttr

      use m_SpatialIntegralV, only : PairedSpatialIntegralsV

      implicit none

! !INPUT PARAMETERS:

      type(AttrVect),    intent(IN) :: inAv1
      type(GeneralGrid), intent(IN) :: GGrid1
      character(len=*),  intent(IN) :: WeightTag1
      type(AttrVect),    intent(IN) :: inAv2
      type(GeneralGrid), intent(IN) :: GGrid2
      character(len=*),  intent(IN) :: WeightTag2
      logical, optional, intent(IN) :: SumWeights
      integer,           intent(IN) :: comm

! !OUTPUT PARAMETERS:

      type(AttrVect),    intent(OUT) :: outAv1
      type(AttrVect),    intent(OUT) :: outAv2

! !REVISION HISTORY:
! 	09May02 - J.W. Larson <larson@mcs.anl.gov> - Initial version.
! 	10Jun02 - J.W. Larson <larson@mcs.anl.gov> - Refactored--now
!                 built on top of PairedIntegralRAttrV_().
!
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::PairedSpatialIntegralRAttrGG_'

       ! Argument Sanity Checks:

  integer :: ierr, length1, length2
  logical :: mySumWeights
  real(FP), dimension(:), pointer :: gridWeights1, gridWeights2

       ! Argument Validity Checks

  if(AttrVect_lsize(inAv1) /= GeneralGrid_lsize(GGrid1)) then
     ierr = AttrVect_lsize(inAv1) - GeneralGrid_lsize(GGrid1)
     write(stderr,'(3a,i8,a,i8)') myname_, &
	  ':: inAv1 / GGrid1 length mismatch:  ', &
	  ' AttrVect_lsize(inAv1) = ',AttrVect_lsize(inAv1), &
	  ' GeneralGrid_lsize(GGrid1) = ',GeneralGrid_lsize(GGrid1)
     call die(myname_)
  endif

  if(AttrVect_lsize(inAv2) /= GeneralGrid_lsize(GGrid2)) then
     ierr = AttrVect_lsize(inAv2) - GeneralGrid_lsize(GGrid2)
     write(stderr,'(3a,i8,a,i8)') myname_, &
	  ':: inAv2 / GGrid2 length mismatch:  ', &
	  ' AttrVect_lsize(inAv2) = ',AttrVect_lsize(inAv2), &
	  ' GeneralGrid_lsize(GGrid2) = ',GeneralGrid_lsize(GGrid2)
     call die(myname_)
  endif

        ! Are we summing the integration weights for either input
        ! GeneralGrid?

  if(present(SumWeights)) then
     mySumWeights = SumWeights
  else
     mySumWeights = .FALSE.
  endif

       ! ensure unambiguous pointer association status for gridWeights1
       ! and gridWeights2

  nullify(gridWeights1) 
  nullify(gridWeights2) 

       ! Extract Grid Weights

  call GeneralGrid_exportRAttr(GGrid1, WeightTag1, gridWeights1, length1)
  call GeneralGrid_exportRAttr(GGrid2, WeightTag2, gridWeights2, length2)


  call PairedSpatialIntegralsV(inAv1, outAv1, gridweights1, WeightTag1, &
	                           inAv2, outAv2, gridweights2, WeightTag2, &
                                   mySumWeights, comm)

       ! Clean up allocated arrays:

  deallocate(gridWeights1, gridWeights2, stat=ierr)
  if(ierr /= 0) then
     write(stderr,'(2a,i8)') myname_, &
	  'ERROR--deallocate(gridWeights1,...) failed, ierr = ',ierr
     call die(myname_)
  endif

 end subroutine PairedSpatialIntegralRAttrGG_

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: PairedSpatialAverageRAttrGG_ - Do two spatial averages at once.
!
! !DESCRIPTION:
! This routine computes spatial averages of the {\tt REAL} attributes
! of the {\tt REAL} attributes of the input {\tt AttrVect} arguments 
! {\tt inAv1} and {\tt inAv2}, returning the integrals in the output 
! {\tt AttrVect} arguments {\tt outAv1} and {\tt outAv2}, respectively .  
! The integrals of {\tt inAv1} and {\tt inAv2} are computed using 
! spatial weights stored in the input {\tt GeneralGrid} arguments  
! {\tt GGrid1} and {\tt GGrid2}, respectively.  The spatial weights in 
! in {\tt GGrid1} and {\tt GGrid2} are identified by the input {\tt CHARACTER} 
! arguments {\tt WeightTag1} and {\tt WeightTag2}, respectively.  
! This paired average is implicitly a 
! distributed operation (the whole motivation for pairing the averages is 
! to reduce communication latency costs), and the Fortran MPI communicator
! handle is defined by the input {\tt INTEGER} argument {\tt comm}.  The 
! summation is an AllReduce operation, with all processes receiving the 
! global sum.
!
! {\bf N.B.:  } The local lengths of the {\tt AttrVect} argument {\tt inAv1} 
! and the {\tt GeneralGrid} {\tt GGrid1} must be equal.  That is, there 
! must be a one-to-one correspondence between the field point values stored 
! in {\tt inAv1} and the point weights stored in {\tt GGrid1}.  The same
! relationship must apply between {\tt inAv2} and {\tt GGrid2}.
!
! {\bf N.B.:  } The output {\tt AttrVect} arguments {\tt outAv1} and 
! {\tt outAv2} are allocated data structures.  The user must deallocate them
!  using the routine {\tt AttrVect\_clean()} when they are no longer needed.  
! Failure to do so will result in a memory leak.
!
! !INTERFACE:

 subroutine PairedSpatialAverageRAttrGG_(inAv1, outAv1, GGrid1, WeightTag1, &
                                         inAv2, outAv2, GGrid2, WeightTag2, &
                                         comm)
! ! USES:

      use m_realkinds, only : FP
   
      use m_stdio
      use m_die
      use m_mpif90

      use m_AttrVect, only : AttrVect
      use m_AttrVect, only : AttrVect_init => init
      use m_AttrVect, only : AttrVect_zero => zero 
      use m_AttrVect, only : AttrVect_clean => clean
      use m_AttrVect, only : AttrVect_lsize => lsize
      use m_AttrVect, only : AttrVect_nRAttr => nRAttr
      use m_AttrVect, only : AttrVect_indexRA => indexRA

      use m_GeneralGrid, only : GeneralGrid
      use m_GeneralGrid, only : GeneralGrid_lsize => lsize
      use m_GeneralGrid, only : GeneralGrid_indexRA => indexRA
      use m_GeneralGrid, only : GeneralGrid_exportRAttr => exportRAttr

      use m_AttrVectReduce, only : AttrVect_LocalWeightedSumRAttr => &
	                                         LocalWeightedSumRAttr

      use m_List, only : List
      use m_List, only : List_nullify => nullify

      implicit none

! !INPUT PARAMETERS:

      type(AttrVect),    intent(IN) :: inAv1
      type(GeneralGrid), intent(IN) :: GGrid1
      character(len=*),  intent(IN) :: WeightTag1
      type(AttrVect),    intent(IN) :: inAv2
      type(GeneralGrid), intent(IN) :: GGrid2
      character(len=*),  intent(IN) :: WeightTag2
      integer,           intent(IN) :: comm

! !OUTPUT PARAMETERS:

      type(AttrVect),    intent(OUT) :: outAv1
      type(AttrVect),    intent(OUT) :: outAv2

! !REVISION HISTORY:
! 	09May02 - J.W. Larson <larson@mcs.anl.gov> - Initial version.
! 	14Jun02 - J.W. Larson <larson@mcs.anl.gov> - Bug fix to reflect
!                 new interface to PairedSpatialIntegralRAttrGG_().
!
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::PairedSpatialAverageRAttrGG_'

  type(AttrVect) :: integratedAv1, integratedAv2
  type(List) :: nullIList
  integer :: i, ierr, iweight1, iweight2

       ! Compute the spatial integral:

     call PairedSpatialIntegralRAttrGG_(inAv1, integratedAv1, GGrid1, WeightTag1, &
	                               inAv2, integratedAv2, GGrid2,     &
				       WeightTag2, .TRUE., comm)


       ! Check value of summed weights (to avoid division by zero):

  iweight1 = AttrVect_indexRA(integratedAv1, WeightTag1)
  if(integratedAv1%rAttr(iweight1, 1) == 0._FP) then   
     write(stderr,'(2a)') myname_, &
	  '::ERROR--Global sum of grid weights in first integral is zero.'
     call die(myname_)
  endif

  iweight2 = AttrVect_indexRA(integratedAv2, WeightTag2)
  if(integratedAv2%rAttr(iweight2, 1) == 0._FP) then
     write(stderr,'(2a)') myname_, &
	  '::ERROR--Global sum of grid weights in second integral is zero.'
     call die(myname_)
  endif

       ! Initialize output AttrVects outAv1 and outAv2:

  call List_nullify(nullIList)

  call AttrVect_init(outAv1, iList=nullIList, rList=inAv1%rList, lsize=1)
  call AttrVect_zero(outAv1)
  call AttrVect_init(outAv2, iList=nullIList, rList=InAv2%rList, lsize=1)
  call AttrVect_zero(outAv2)

       ! Divide by global weight sum to compute spatial averages from 
       ! spatial integrals.

  do i=1,AttrVect_nRAttr(outAv1)
     outAv1%rAttr(i,1) = integratedAv1%rAttr(i,1) &
	                               / integratedAv1%rAttr(iweight1,1) 
  end do

  do i=1,AttrVect_nRAttr(outAv2)
     outAv2%rAttr(i,1) = integratedAv2%rAttr(i,1) &
	                               / integratedAv2%rAttr(iweight2,1) 
  end do

       ! Clean up temporary AttrVects:

  call AttrVect_clean(integratedAv1)
  call AttrVect_clean(integratedAv2)

 end subroutine PairedSpatialAverageRAttrGG_


!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: PairedMaskedIntegralRAttrGG_ - Do two masked integrals at once.
!
! !DESCRIPTION:
! This routine computes a pair of masked spatial integrals of the {\tt REAL} 
! attributes of the input {\tt AttrVect} arguments {\tt inAv} and 
! {\tt inAv2}, returning the masked integrals in the output {\tt AttrVect} 
! {\tt outAv1} and {\tt outAv2}, respectively.  All of the spatial weighting
! and masking data for each set of integrals are assumed stored in the input 
! {\tt GeneralGrid} arguments {\tt GGrid} and {\tt GGrid2}.  If integer 
! masks are to be used, their integer  attribute names in {\tt GGrid1} 
! and {\tt GGrid2} are named as a colon-delimited lists in the optional 
! {\tt CHARACTER} input arguments {\tt iMaskTags1} and {\tt iMaskTags2}, 
! respectively.  Real masks (if desired) are referenced by their real 
! attribute names in {\tt GGrid1} and {\tt GGrid2} are named as 
! colon-delimited lists in the optional {\tt CHARACTER} input arguments 
! {\tt rMaskTags1} and {\tt rMaskTags2}, respectively.  The user specifies 
! a choice of mask combination method with the input {\tt LOGICAL} argument
! {\tt UseFastMethod}.  If ${\tt UseFastMethod} = {\tt .FALSE.}$ this 
! routine checks each mask entry to ensure that the integer masks contain
! only ones and zeroes, and that entries in the real masks are all in
! the closed interval $[0,1]$.  If ${\tt UseFastMethod} = {\tt .TRUE.}$,
! this routine performs direct products of the masks, assuming that the
! user has validated them in advance.  The optional {\tt LOGICAL} input 
! argument {\tt SumWeights} determines whether the masked sum of the spatial
! weights is computed and returned in {\tt outAv1} and {\tt outAv2} with the 
! real attribute names supplied in the {\tt CHARACTER} input arguments
! {\tt SpatialWeightTag1}, and {\tt SpatialWeightTag2}, respectively.  
! This paired integral is implicitly a distributed operation (the whole 
! motivation for pairing the averages is to reduce communication latency 
! costs), and the Fortran MPI communicator handle is defined by the input 
! {\tt INTEGER} argument {\tt comm}.  The 
! summation is an AllReduce operation, with all processes receiving the 
! global sum.
!
! {\bf N.B.:  } The local lengths of the {\tt AttrVect} argument {\tt inAv1} 
! and the {\tt GeneralGrid} {\tt GGrid1} must be equal.  That is, there 
! must be a one-to-one correspondence between the field point values stored 
! in {\tt inAv1} and the point weights stored in {\tt GGrid1}.  The same
! relationship must apply between {\tt inAv2} and {\tt GGrid2}.
!
! {\bf N.B.:  }  If {\tt PairedMaskedIntegralRAttrGG\_()} is invoked with the 
! optional {\tt LOGICAL} input argument {\tt SumWeights} set as {\tt .TRUE.}, 
! then the value of {\tt SpatialWeightTag1} must not conflict with any of the 
! {\tt REAL} attribute tags in {\tt inAv1} and the value of 
! {\tt SpatialWeightTag2} must not conflict with any of the {\tt REAL} 
! attribute tags in {\tt inAv2}.
!
! {\bf N.B.:  } The output {\tt AttrVect} arguments {\tt outAv1} and 
! {\tt outAv2} are allocated data structures.  The user must deallocate them
!  using the routine {\tt AttrVect\_clean()} when they are no longer needed.  
! Failure to do so will result in a memory leak.
!
! !INTERFACE:

 subroutine PairedMaskedIntegralRAttrGG_(inAv1, outAv1, GGrid1, &
                                         SpatialWeightTag1, rMaskTags1, &
                                         iMaskTags1, inAv2, outAv2, GGrid2, &
                                         SpatialWeightTag2, rMaskTags2, &
                                         iMaskTags2, UseFastMethod, &
                                         SumWeights, comm)
! ! USES:

      use m_stdio
      use m_die
      use m_mpif90

      use m_realkinds, only : FP

      use m_AttrVect, only : AttrVect
      use m_AttrVect, only : AttrVect_lsize => lsize
      use m_AttrVect, only : AttrVect_nRAttr => nRAttr

      use m_GeneralGrid, only : GeneralGrid
      use m_GeneralGrid, only : GeneralGrid_lsize => lsize
      use m_GeneralGrid, only : GeneralGrid_indexRA => indexRA
      use m_GeneralGrid, only : GeneralGrid_exportRAttr => exportRAttr

      use m_AttrVectReduce, only : AttrVect_LocalWeightedSumRAttr => &
	                                         LocalWeightedSumRAttr

      implicit none

! !INPUT PARAMETERS:

      type(AttrVect),              intent(IN) :: inAv1
      type(GeneralGrid),           intent(IN) :: GGrid1
      character(len=*),            intent(IN) :: SpatialWeightTag1
      character(len=*),  optional, intent(IN) :: iMaskTags1
      character(len=*),  optional, intent(IN) :: rMaskTags1
      type(AttrVect),              intent(IN) :: inAv2
      type(GeneralGrid),           intent(IN) :: GGrid2
      character(len=*),            intent(IN) :: SpatialWeightTag2
      character(len=*),  optional, intent(IN) :: iMaskTags2
      character(len=*),  optional, intent(IN) :: rMaskTags2
      logical,                     intent(IN) :: UseFastMethod
      logical,           optional, intent(IN) :: SumWeights
      integer,                     intent(IN) :: comm

! !OUTPUT PARAMETERS:

      type(AttrVect),    intent(OUT) :: outAv1
      type(AttrVect),    intent(OUT) :: outAv2

! !REVISION HISTORY:
! 	17Jun02 - J.W. Larson <larson@mcs.anl.gov> - Initial version.
! 	19Jun02 - J.W. Larson <larson@mcs.anl.gov> - Shortened the name
!                 for compatibility with the Portland Group f90 compiler
!
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_ = &
                            myname//'::PairedMaskedIntegralRAttrGG_'

  logical :: mySumWeights
  real(FP), dimension(:), pointer :: PairedBuffer, OutPairedBuffer
  integer :: ierr, nRA1, nRA2, PairedBufferLength

        ! Basic Argument Validity Checks:

  if(AttrVect_lsize(inAv1) /= GeneralGrid_lsize(GGrid1)) then
     ierr = AttrVect_lsize(inAv1) - GeneralGrid_lsize(GGrid1)
     write(stderr,'(3a,i8,a,i8)') myname_, &
          ':: inAv1 / GGrid1 length mismatch:  ', &
          ' AttrVect_lsize(inAv1) = ',AttrVect_lsize(inAv1), &
          ' GeneralGrid_lsize(GGrid1) = ',GeneralGrid_lsize(GGrid1)
     call die(myname_)
  endif

  if(AttrVect_lsize(inAv2) /= GeneralGrid_lsize(GGrid2)) then
     ierr = AttrVect_lsize(inAv2) - GeneralGrid_lsize(GGrid2)
     write(stderr,'(3a,i8,a,i8)') myname_, &
          ':: inAv2 / GGrid2 length mismatch:  ', &
          ' AttrVect_lsize(inAv2) = ',AttrVect_lsize(inAv2), &
          ' GeneralGrid_lsize(GGrid2) = ',GeneralGrid_lsize(GGrid2)
     call die(myname_)
  endif

        ! Are we summing the integration weights for the input
        ! GeneralGrids?

  if(present(SumWeights)) then
     mySumWeights = SumWeights
  else
     mySumWeights = .FALSE.
  endif

        ! Begin by invoking MaskedSpatialIntegralRAttrGG_() for each
        ! AttrVect/GeneralGrid pair.  This is done LOCALLY to create
        ! integratedAv1 and integratedAv2, respectively.

        ! Local Masked Integral #1:

  if(present(iMaskTags1)) then

     if(present(rMaskTags1)) then ! both Integer and Real Masking
	call MaskedSpatialIntegralRAttrGG_(inAv1, outAv1, GGrid1, &
	                                  SpatialWeightTag1, iMaskTags1, &
					  rMaskTags1, UseFastMethod,  &
					  mySumWeights, SpatialWeightTag1)
     else ! Integer Masking Only
	call MaskedSpatialIntegralRAttrGG_(inAv1, outAv1, GGrid1, &
	                                  SpatialWeightTag1, &
					  iMaskTags=iMaskTags1, &
					  UseFastMethod=UseFastMethod,  &
					  SumWeights=mySumWeights, &
					  WeightSumTag=SpatialWeightTag1)
     endif ! if(present(rMaskTags1))...

  else ! No Integer Masking

     if(present(rMaskTags1)) then ! Real Masking Only
	call MaskedSpatialIntegralRAttrGG_(inAv1, outAv1, GGrid1, &
	                                  SpatialWeightTag=SpatialWeightTag1, &
					  rMaskTags=rMaskTags1, &
					  UseFastMethod=UseFastMethod,  &
					  SumWeights=mySumWeights, &
					  WeightSumTag=SpatialWeightTag1)
     else ! Neither Integer nor Real Masking
	call MaskedSpatialIntegralRAttrGG_(inAv1, outAv1, GGrid1, &
	                                  SpatialWeightTag=SpatialWeightTag1, &
					  UseFastMethod=UseFastMethod,  &
					  SumWeights=mySumWeights, &
					  WeightSumTag=SpatialWeightTag1)

     endif ! if(present(rMaskTags1))...

  endif ! if(present(iMaskTags1))...

        ! Local Masked Integral #2:

  if(present(iMaskTags2)) then

     if(present(rMaskTags2)) then ! both Integer and Real Masking
	call MaskedSpatialIntegralRAttrGG_(inAv2, outAv2, GGrid2, &
	                                  SpatialWeightTag2, iMaskTags2, &
					  rMaskTags2, UseFastMethod,  &
					  mySumWeights, SpatialWeightTag2)
     else ! Integer Masking Only
	call MaskedSpatialIntegralRAttrGG_(inAv2, outAv2, GGrid2, &
	                                  SpatialWeightTag2, &
					  iMaskTags=iMaskTags2, &
					  UseFastMethod=UseFastMethod,  &
					  SumWeights=mySumWeights, &
					  WeightSumTag=SpatialWeightTag2)
     endif ! if(present(rMaskTags2))...

  else ! No Integer Masking

     if(present(rMaskTags2)) then ! Real Masking Only
	call MaskedSpatialIntegralRAttrGG_(inAv2, outAv2, GGrid2, &
	                                  SpatialWeightTag=SpatialWeightTag2, &
					  rMaskTags=rMaskTags2, &
					  UseFastMethod=UseFastMethod,  &
					  SumWeights=mySumWeights, &
					  WeightSumTag=SpatialWeightTag2)
     else ! Neither Integer nor Real Masking
	call MaskedSpatialIntegralRAttrGG_(inAv2, outAv2, GGrid2, &
	                                  SpatialWeightTag=SpatialWeightTag2, &
					  UseFastMethod=UseFastMethod,  &
					  SumWeights=mySumWeights, &
					  WeightSumTag=SpatialWeightTag2)

     endif ! if(present(rMaskTags2))...

  endif ! if(present(iMaskTags2))...

       ! Create the paired buffer for the Global Sum

  nRA1 = AttrVect_nRAttr(outAv1)
  nRA2 = AttrVect_nRAttr(outAv2)

  PairedBufferLength =  nRA1 + nRA2
  allocate(PairedBuffer(PairedBufferLength), OutPairedBuffer(PairedBufferLength), &
           stat=ierr)
  if(ierr /= 0) then
     write(stderr,'(2a,i8)') myname_, &
	  ':: Fatal error--allocate(PairedBuffer...failed, ierr = ',ierr
     call die(myname_)
  endif

       ! Load the paired buffer

  PairedBuffer(1:nRA1) = outAv1%rAttr(1:nRA1,1)
  PairedBuffer(nRA1+1:PairedBufferLength) = outAv2%rAttr(1:nRA2,1)

       ! Perform the global sum on the paired buffer

  call MPI_AllReduce(PairedBuffer, OutPairedBuffer, PairedBufferLength, &
                        MP_Type(PairedBuffer(1)), MP_SUM, comm, ierr)
  if(ierr /= 0) then
     write(stderr,'(2a,i8)') myname_, &
  	      ':: Fatal Error--MPI_ALLREDUCE() failed with ierror = ',ierr
     call MP_perr_die(myname_,'MPI_ALLREDUCE() failed',ierr)
  endif

       ! Unload OutPairedBuffer into outAv1 and outAv2:

  outAv1%rAttr(1:nRA1,1) = OutPairedBuffer(1:nRA1)
  outAv2%rAttr(1:nRA2,1) = OutPairedBuffer(nRA1+1:PairedBufferLength)

  deallocate(PairedBuffer, OutPairedBuffer, stat=ierr)
  if(ierr /= 0) then
     write(stderr,'(2a,i8)') myname_, &
	  ':: Fatal error--deallocate(PairedBuffer...failed, ierr = ',ierr
     call die(myname_)
  endif

 end subroutine PairedMaskedIntegralRAttrGG_

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: PairedMaskedAverageRAttrGG_ - Do two masked averages at once.
!
! !DESCRIPTION:
! This routine computes a pair of masked spatial averages of the {\tt REAL} 
! attributes of the input {\tt AttrVect} arguments {\tt inAv} and 
! {\tt inAv2}, returning the masked averagess in the output {\tt AttrVect} 
! {\tt outAv1} and {\tt outAv2}, respectively.  All of the spatial weighting
! and masking data for each set of averages are assumed stored in the input 
! {\tt GeneralGrid} arguments {\tt GGrid} and {\tt GGrid2}.  If integer 
! masks are to be used, their integer  attribute names in {\tt GGrid1} 
! and {\tt GGrid2} are named as a colon-delimited lists in the optional 
! {\tt CHARACTER} input arguments {\tt iMaskTags1} and {\tt iMaskTags2}, 
! respectively.  Real masks (if desired) are referenced by their real 
! attribute names in {\tt GGrid1} and {\tt GGrid2} are named as 
! colon-delimited lists in the optional {\tt CHARACTER} input arguments 
! {\tt rMaskTags1} and {\tt rMaskTags2}, respectively.  The user specifies 
! a choice of mask combination method with the input {\tt LOGICAL} argument
! {\tt UseFastMethod}.  If ${\tt UseFastMethod} = {\tt .FALSE.}$ this 
! routine checks each mask entry to ensure that the integer masks contain
! only ones and zeroes, and that entries in the real masks are all in
! the closed interval $[0,1]$.  If ${\tt UseFastMethod} = {\tt .TRUE.}$,
! this routine performs direct products of the masks, assuming that the
! user has validated them in advance.  This paired average is implicitly 
! a distributed operation (the whole motivation for pairing the averages 
! is to reduce communication latency costs), and the Fortran MPI communicator 
! handle is defined by the input {\tt INTEGER} argument {\tt comm}.  The 
! summation is an AllReduce operation, with all processes receiving the 
! global sum.
!
! {\bf N.B.:  } The local lengths of the {\tt AttrVect} argument {\tt inAv1} 
! and the {\tt GeneralGrid} {\tt GGrid1} must be equal.  That is, there 
! must be a one-to-one correspondence between the field point values stored 
! in {\tt inAv1} and the point weights stored in {\tt GGrid1}.  The same
! relationship must apply between {\tt inAv2} and {\tt GGrid2}.
!
! {\bf N.B.:  } The output {\tt AttrVect} arguments {\tt outAv1} and 
! {\tt outAv2} are allocated data structures.  The user must deallocate them
!  using the routine {\tt AttrVect\_clean()} when they are no longer needed.  
! Failure to do so will result in a memory leak.
!
! !INTERFACE:

 subroutine PairedMaskedAverageRAttrGG_(inAv1, outAv1, GGrid1, &
                                        SpatialWeightTag1, rMaskTags1, &
                                        iMaskTags1, inAv2, outAv2, GGrid2, &
                                        SpatialWeightTag2, rMaskTags2, &
                                        iMaskTags2, UseFastMethod, &
                                        comm)
! ! USES:

      use m_stdio
      use m_die
      use m_mpif90

      use m_realkinds, only : FP

      use m_AttrVect, only : AttrVect
      use m_AttrVect, only : AttrVect_init => init
      use m_AttrVect, only : AttrVect_zero => zero 
      use m_AttrVect, only : AttrVect_clean => clean
      use m_AttrVect, only : AttrVect_lsize => lsize
      use m_AttrVect, only : AttrVect_nRAttr => nRAttr

      use m_GeneralGrid, only : GeneralGrid
      use m_GeneralGrid, only : GeneralGrid_lsize => lsize
      use m_GeneralGrid, only : GeneralGrid_indexRA => indexRA
      use m_GeneralGrid, only : GeneralGrid_exportRAttr => exportRAttr

      use m_AttrVectReduce, only : AttrVect_LocalWeightedSumRAttr => &
	                                         LocalWeightedSumRAttr

      use m_List, only : List
      use m_List, only : List_nullify => nullify

      implicit none

! !INPUT PARAMETERS:

      type(AttrVect),              intent(IN) :: inAv1
      type(GeneralGrid),           intent(IN) :: GGrid1
      character(len=*),            intent(IN) :: SpatialWeightTag1
      character(len=*),  optional, intent(IN) :: iMaskTags1
      character(len=*),  optional, intent(IN) :: rMaskTags1
      type(AttrVect),              intent(IN) :: inAv2
      type(GeneralGrid),           intent(IN) :: GGrid2
      character(len=*),            intent(IN) :: SpatialWeightTag2
      character(len=*),  optional, intent(IN) :: iMaskTags2
      character(len=*),  optional, intent(IN) :: rMaskTags2
      logical,                     intent(IN) :: UseFastMethod
      integer,                     intent(IN) :: comm

! !OUTPUT PARAMETERS:

      type(AttrVect),    intent(OUT) :: outAv1
      type(AttrVect),    intent(OUT) :: outAv2

! !REVISION HISTORY:
! 	17Jun02 - J.W. Larson <larson@mcs.anl.gov> - Initial version.
! 	19Jun02 - J.W. Larson <larson@mcs.anl.gov> - Shortened the name
!                 for compatibility with the Portland Group f90 compiler
! 	25Jul02 - J.W. Larson E.T. Ong - Bug fix.  This routine was 
!                 previously doing integrals rather than area averages.
!
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_ = &
                            myname//'::PairedMaskedAverageRAttrGG_'

  type(AttrVect) :: LocalIntegral1, LocalIntegral2
  type(List) :: nullIList
  real(FP), dimension(:), pointer :: PairedBuffer, OutPairedBuffer
  integer :: i, ierr, nRA1, nRA2, PairedBufferLength
  real(FP) :: WeightSumInv

        ! Basic Argument Validity Checks:

  if(AttrVect_lsize(inAv1) /= GeneralGrid_lsize(GGrid1)) then
     ierr = AttrVect_lsize(inAv1) - GeneralGrid_lsize(GGrid1)
     write(stderr,'(3a,i8,a,i8)') myname_, &
          ':: inAv1 / GGrid1 length mismatch:  ', &
          ' AttrVect_lsize(inAv1) = ',AttrVect_lsize(inAv1), &
          ' GeneralGrid_lsize(GGrid1) = ',GeneralGrid_lsize(GGrid1)
     call die(myname_)
  endif

  if(AttrVect_lsize(inAv2) /= GeneralGrid_lsize(GGrid2)) then
     ierr = AttrVect_lsize(inAv2) - GeneralGrid_lsize(GGrid2)
     write(stderr,'(3a,i8,a,i8)') myname_, &
          ':: inAv2 / GGrid2 length mismatch:  ', &
          ' AttrVect_lsize(inAv2) = ',AttrVect_lsize(inAv2), &
          ' GeneralGrid_lsize(GGrid2) = ',GeneralGrid_lsize(GGrid2)
     call die(myname_)
  endif

        ! Begin by invoking MaskedSpatialIntegralRAttrGG_() for each
        ! AttrVect/GeneralGrid pair.  This is done LOCALLY to create
        ! LocalIntegral1 and LocalIntegral2, respectively.

        ! Local Masked Integral #1:

  if(present(iMaskTags1)) then

     if(present(rMaskTags1)) then ! both Integer and Real Masking
	call MaskedSpatialIntegralRAttrGG_(inAv1, LocalIntegral1, GGrid1, &
	                                  SpatialWeightTag1, iMaskTags1, &
					  rMaskTags1, UseFastMethod,  &
					  .TRUE., SpatialWeightTag1)
     else ! Integer Masking Only
	call MaskedSpatialIntegralRAttrGG_(inAv1, LocalIntegral1, GGrid1, &
	                                  SpatialWeightTag1, &
					  iMaskTags=iMaskTags1, &
					  UseFastMethod=UseFastMethod,  &
					  SumWeights=.TRUE., &
					  WeightSumTag=SpatialWeightTag1)
     endif ! if(present(rMaskTags1))...

  else ! No Integer Masking

     if(present(rMaskTags1)) then ! Real Masking Only
	call MaskedSpatialIntegralRAttrGG_(inAv1, LocalIntegral1, GGrid1, &
	                                  SpatialWeightTag=SpatialWeightTag1, &
					  rMaskTags=rMaskTags1, &
					  UseFastMethod=UseFastMethod,  &
					  SumWeights=.TRUE., &
					  WeightSumTag=SpatialWeightTag1)
     else ! Neither Integer nor Real Masking
	call MaskedSpatialIntegralRAttrGG_(inAv1, LocalIntegral1, GGrid1, &
	                                  SpatialWeightTag=SpatialWeightTag1, &
					  UseFastMethod=UseFastMethod,  &
					  SumWeights=.TRUE., &
					  WeightSumTag=SpatialWeightTag1)

     endif ! if(present(rMaskTags1))...

  endif ! if(present(iMaskTags1))...

        ! Local Masked Integral #2:

  if(present(iMaskTags2)) then

     if(present(rMaskTags2)) then ! both Integer and Real Masking
	call MaskedSpatialIntegralRAttrGG_(inAv2, LocalIntegral2, GGrid2, &
	                                  SpatialWeightTag2, iMaskTags2, &
					  rMaskTags2, UseFastMethod,  &
					  .TRUE., SpatialWeightTag2)
     else ! Integer Masking Only
	call MaskedSpatialIntegralRAttrGG_(inAv2, LocalIntegral2, GGrid2, &
	                                  SpatialWeightTag2, &
					  iMaskTags=iMaskTags2, &
					  UseFastMethod=UseFastMethod,  &
					  SumWeights=.TRUE., &
					  WeightSumTag=SpatialWeightTag2)
     endif ! if(present(rMaskTags2))...

  else ! No Integer Masking

     if(present(rMaskTags2)) then ! Real Masking Only
	call MaskedSpatialIntegralRAttrGG_(inAv2, LocalIntegral2, GGrid2, &
	                                  SpatialWeightTag=SpatialWeightTag2, &
					  rMaskTags=rMaskTags2, &
					  UseFastMethod=UseFastMethod,  &
					  SumWeights=.TRUE., &
					  WeightSumTag=SpatialWeightTag2)
     else ! Neither Integer nor Real Masking
	call MaskedSpatialIntegralRAttrGG_(inAv2, LocalIntegral2, GGrid2, &
	                                  SpatialWeightTag=SpatialWeightTag2, &
					  UseFastMethod=UseFastMethod,  &
					  SumWeights=.TRUE., &
					  WeightSumTag=SpatialWeightTag2)

     endif ! if(present(rMaskTags2))...

  endif ! if(present(iMaskTags2))...

       ! Create the paired buffer for the Global Sum

  nRA1 = AttrVect_nRAttr(LocalIntegral1)
  nRA2 = AttrVect_nRAttr(LocalIntegral2)

  PairedBufferLength =  nRA1 + nRA2
  allocate(PairedBuffer(PairedBufferLength), OutPairedBuffer(PairedBufferLength), &
           stat=ierr)
  if(ierr /= 0) then
     write(stderr,'(2a,i8)') myname_, &
	  ':: Fatal error--allocate(PairedBuffer...failed, ierr = ',ierr
     call die(myname_)
  endif

       ! Load the paired buffer

  PairedBuffer(1:nRA1) = LocalIntegral1%rAttr(1:nRA1,1)
  PairedBuffer(nRA1+1:PairedBufferLength) = LocalIntegral2%rAttr(1:nRA2,1)

       ! Perform the global sum on the paired buffer

  call MPI_AllReduce(PairedBuffer, OutPairedBuffer, PairedBufferLength, &
                        MP_Type(PairedBuffer(1)), MP_SUM, comm, ierr)
  if(ierr /= 0) then
     write(stderr,'(2a,i8)') myname_, &
  	      ':: Fatal Error--MPI_ALLREDUCE() failed with ierror = ',ierr
     call MP_perr_die(myname_,'MPI_ALLREDUCE() failed',ierr)
  endif

       ! Create outAv1 and outAv2 from inAv1 and inAv2, respectively:

  call List_nullify(nullIList)

  call AttrVect_init(outAv1, iList=nullIList, rList=inAv1%rList, lsize=1)
  call AttrVect_zero(outAv1)
  call AttrVect_init(outAv2, iList=nullIList, rList=inAv2%rList, lsize=1)
  call AttrVect_zero(outAv2)

       ! Unload/rescale OutPairedBuffer into outAv1 and outAv2:

  nRA1 = AttrVect_nRAttr(outAv1)
  nRA2 = AttrVect_nRAttr(outAv2)

       ! First outAv1:

  if(OutPairedBuffer(nRA1+1) /= 0.) then
     WeightSumInv = 1._FP / OutPairedBuffer(nRA1+1) ! Sum of weights on grid1
                                                 ! is the nRA1+1th element in
                                                 ! the paired buffer.
  else
     write(stderr,'(2a)') myname_, &
	  ':: FATAL ERROR--Sum of the Weights for integral #1 is zero!  Terminating...'
     call die(myname_)
  endif

       ! Rescale global integral to get global average:

  do i=1,nRA1
     outAv1%rAttr(i,1) = WeightSumInv * OutPairedBuffer(i)
  end do

       ! And then outAv2:

  if(OutPairedBuffer(PairedBufferLength) /= 0.) then
     WeightSumInv = 1._FP / OutPairedBuffer(PairedBufferLength) ! Sum of weights on grid2
                                                             ! is the last element in
                                                             ! the paired buffer.
  else
     write(stderr,'(2a)') myname_, &
	  ':: FATAL ERROR--Sum of the Weights for integral #2 is zero!  Terminating...'
     call die(myname_)
  endif

       ! Rescale global integral to get global average:

  do i=1,nRA2
     outAv2%rAttr(i,1) = WeightSumInv * OutPairedBuffer(i+nRA1+1)
  end do

       ! Clean up allocated structures

  call AttrVect_clean(LocalIntegral1)
  call AttrVect_clean(LocalIntegral2)

  deallocate(PairedBuffer, OutPairedBuffer, stat=ierr)
  if(ierr /= 0) then
     write(stderr,'(2a,i8)') myname_, &
	  ':: Fatal error--deallocate(PairedBuffer...failed, ierr = ',ierr
     call die(myname_)
  endif

 end subroutine PairedMaskedAverageRAttrGG_

 end module m_SpatialIntegral