Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

m_SpatialIntegralV.F90 @ 4775

Last change on this file since 4775 was 4775, checked in by aclsce, 5 years ago

Imported oasis3-mct from Cerfacs svn server (not suppotred anymore).

The version has been extracted from https://oasis3mct.cerfacs.fr/svn/branches/OASIS3-MCT_2.0_branch/oasis3-mct@1818

File size: 69.3 KB

Line
1	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2	! Math and Computer Science Division, Argonne National Laboratory !
3	!-----------------------------------------------------------------------
4	! CVS m_SpatialIntegralV.F90,v 1.7 2008-05-12 01:59:33 jacob Exp
5	! CVS MCT_2_8_0
6	!BOP -------------------------------------------------------------------
7	!
8	! !MODULE: m_SpatialIntegralV - Spatial Integrals and Averages using vectors of weights
9	!
10	! !DESCRIPTION: This module provides spatial integration and averaging
11	! services for the MCT similar to those in {\tt m\_SpatialIntegral} except
12	! the weights are provided by an input vector instead of through a
13	! {\tt GeneralGrid}. See the description for {\tt m\_SpatialIntegral} for
14	! more information
15	!
16	!
17	! Paired masked spatial integrals and averages have not yet been implemented in
18	! vector form.
19	!
20	! !INTERFACE:
21
22	module m_SpatialIntegralV
23
24	implicit none
25
26	private ! except
27
28	! !PUBLIC MEMBER FUNCTIONS:
29
30	public :: SpatialIntegralV ! Spatial Integral
31	public :: SpatialAverageV ! Spatial Area Average
32
33	public :: MaskedSpatialIntegralV ! Masked Spatial Integral
34	public :: MaskedSpatialAverageV ! MaskedSpatial Area Average
35
36	public :: PairedSpatialIntegralsV ! A Pair of Spatial
37	! Integrals
38
39	public :: PairedSpatialAveragesV ! A Pair of Spatial
40	! Area Averages
41
42	interface SpatialIntegralV ; module procedure &
43	SpatialIntegralRAttrVSP_, &
44	SpatialIntegralRAttrVDP_
45	end interface
46	interface SpatialAverageV ; module procedure &
47	SpatialAverageRAttrVSP_, &
48	SpatialAverageRAttrVDP_
49	end interface
50	interface MaskedSpatialIntegralV ; module procedure &
51	MaskedSpatialIntegralRAttrVSP_, &
52	MaskedSpatialIntegralRAttrVDP_
53	end interface
54	interface MaskedSpatialAverageV ; module procedure &
55	MaskedSpatialAverageRAttrVSP_, &
56	MaskedSpatialAverageRAttrVDP_
57	end interface
58	interface PairedSpatialIntegralsV ; module procedure &
59	PairedSpatialIntegralRAttrVSP_, &
60	PairedSpatialIntegralRAttrVDP_
61	end interface
62	interface PairedSpatialAveragesV ; module procedure &
63	PairedSpatialAverageRAttrVSP_, &
64	PairedSpatialAverageRAttrVDP_
65	end interface
66
67	! !REVISION HISTORY:
68	! 4Jan04 - R.Jacob <jacob@mcs.anl.gov> - move Vector versions of routines
69	! from m_SpatialIntegral to this file.
70	!EOP ___________________________________________________________________
71
72	character(len=*),parameter :: myname='MCT::m_SpatialIntegralV'
73
74	contains
75
76	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
77	! Math and Computer Science Division, Argonne National Laboratory !
78	!BOP -------------------------------------------------------------------
79	!
80	! !IROUTINE: SpatialIntegralRAttrVSP_ - Compute spatial integral.
81	!
82	! !DESCRIPTION:
83	! This routine computes spatial integrals of the {\tt REAL} attributes
84	! of the {\tt REAL} attributes of the input {\tt AttrVect} argument
85	! {\tt inAv}. {\tt SpatialIntegralRAttrV\_()} takes the input
86	! {\tt AttrVect} argument {\tt inAv} and computes the spatial
87	! integral using weights stored in the input {\tt REAL} array argument
88	! {\tt Weights}. The integral of each {\tt REAL} attribute is returned
89	! in the output {\tt AttrVect} argument {\tt outAv}. If
90	! {\tt SpatialIntegralRAttrV\_()} is invoked with the optional {\tt LOGICAL}
91	! input argument {\tt SumWeights} set as {\tt .TRUE.}, then the weights
92	! are also summed and stored in {\tt outAv} (and can be referenced with
93	! the attribute name {\tt WeightTag}. If {\tt SpatialIntegralRAttrV\_()} is
94	! invoked with the optional {\tt INTEGER} argument {\tt comm} (a Fortran
95	! MPI communicator handle), the summation operations for the integral are
96	! completed on the local process, then reduced across the communicator,
97	! with all processes receiving the result.
98	!
99	! {\bf N.B.: } The local lengths of the {\tt AttrVect} argument {\tt inAv}
100	! and the input array {\tt Weights} must be equal. That is, there must be
101	! a one-to-one correspondence between the field point values stored
102	! in {\tt inAv} and the point weights stored in {\tt Weights}.
103	!
104	! {\bf N.B.: } If {\tt SpatialIntegralRAttrV\_()} is invoked with the
105	! optional {\tt LOGICAL} input argument {\tt SumWeights} set as {\tt .TRUE.}.
106	! In this case, the none of {\tt REAL} attribute tags in {\tt inAv} may be
107	! named the same as the string contained in {\tt WeightTag}, which is an
108	! attribute name reserved for the sum of the weights in the output {\tt AttrVect}
109	! {\tt outAv}.
110	!
111	! {\bf N.B.: } The output {\tt AttrVect} argument {\tt outAv} is an
112	! allocated data structure. The user must deallocate it using the routine
113	! {\tt AttrVect\_clean()} when it is no longer needed. Failure to do so
114	! will result in a memory leak.
115	!
116	! !INTERFACE:
117
118	subroutine SpatialIntegralRAttrVSP_(inAv, outAv, Weights, SumWeights, &
119	WeightTag, comm)
120
121	! ! USES:
122
123	use m_stdio
124	use m_die
125	use m_mpif90
126	use m_realkinds, only : SP
127
128	use m_AttrVect, only : AttrVect
129	use m_AttrVect, only : AttrVect_lsize => lsize
130
131	use m_AttrVectReduce, only : AttrVect_GlobalWeightedSumRAttr => &
132	GlobalWeightedSumRAttr
133	use m_AttrVectReduce, only : AttrVect_LocalWeightedSumRAttr => &
134	LocalWeightedSumRAttr
135
136	implicit none
137
138	! !INPUT PARAMETERS:
139
140	type(AttrVect), intent(IN) :: inAv
141	real(SP), dimension(:), pointer :: Weights
142	logical, optional, intent(IN) :: SumWeights
143	character(len=*), optional, intent(IN) :: WeightTag
144	integer, optional, intent(IN) :: comm
145
146	! !OUTPUT PARAMETERS:
147
148	type(AttrVect), intent(OUT) :: outAv
149
150	! !REVISION HISTORY:
151	! 07Jun02 - J.W. Larson <larson@mcs.anl.gov> - initial version
152	!EOP ___________________________________________________________________
153
154	character(len=*),parameter :: myname_=myname//'::SpatialIntegralRAttrVSP_'
155
156	integer :: ierr, length
157	logical :: mySumWeights
158
159	! Argument Validity Checks
160
161	if(AttrVect_lsize(inAv) /= size(Weights)) then
162	ierr = AttrVect_lsize(inAv) - size(Weights)
163	write(stderr,'(3a,i8,a,i8)') myname_, &
164	':: inAv / Weights array length mismatch: ', &
165	' AttrVect_lsize(inAv) = ',AttrVect_lsize(inAv), &
166	' size(Weights) = ',size(Weights)
167	call die(myname_)
168	endif
169
170	if(present(SumWeights)) then
171	mySumWeights = SumWeights
172	if(.not. present(WeightTag)) then
173	write(stderr,'(3a)') myname_,':: FATAL--If the input argument SumWeights=.TRUE.,', &
174	' then the argument WeightTag must be provided.'
175	call die(myname_)
176	endif
177	else
178	mySumWeights = .FALSE.
179	endif
180
181	! Compute the sum
182
183	if(present(comm)) then ! compute distributed AllReduce-style sum:
184
185	if(mySumWeights) then ! return the spatial sum of the weights in outAV
186	call AttrVect_GlobalWeightedSumRAttr(inAV, outAV, Weights, &
187	comm, WeightTag)
188	else
189	call AttrVect_GlobalWeightedSumRAttr(inAV, outAV, Weights, comm)
190	endif
191
192	else ! compute local sum:
193
194	if(mySumWeights) then ! return the spatial sum of the weights in outAV
195	call AttrVect_LocalWeightedSumRAttr(inAV, outAV, Weights, &
196	WeightTag)
197	else
198	call AttrVect_LocalWeightedSumRAttr(inAV, outAV, Weights)
199	endif
200
201	endif ! if(present(comm))...
202
203	end subroutine SpatialIntegralRAttrVSP_
204
205	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
206	! Math and Computer Science Division, Argonne National Laboratory !
207	! -------------------------------------------------------------------
208	!
209	! !IROUTINE: SpatialIntegralRAttrVDP_ - Compute spatial integral.
210	!
211	! !DESCRIPTION:
212	! Double precision version of SpatialIntegralRAttrVSP_
213	!
214	! !INTERFACE:
215
216	subroutine SpatialIntegralRAttrVDP_(inAv, outAv, Weights, SumWeights, &
217	WeightTag, comm)
218
219	! ! USES:
220
221	use m_stdio
222	use m_die
223	use m_mpif90
224	use m_realkinds, only : DP
225
226	use m_AttrVect, only : AttrVect
227	use m_AttrVect, only : AttrVect_lsize => lsize
228
229	use m_AttrVectReduce, only : AttrVect_GlobalWeightedSumRAttr => &
230	GlobalWeightedSumRAttr
231	use m_AttrVectReduce, only : AttrVect_LocalWeightedSumRAttr => &
232	LocalWeightedSumRAttr
233
234	implicit none
235
236	! !INPUT PARAMETERS:
237
238	type(AttrVect), intent(IN) :: inAv
239	real(DP), dimension(:), pointer :: Weights
240	logical, optional, intent(IN) :: SumWeights
241	character(len=*), optional, intent(IN) :: WeightTag
242	integer, optional, intent(IN) :: comm
243
244	! !OUTPUT PARAMETERS:
245
246	type(AttrVect), intent(OUT) :: outAv
247
248	! !REVISION HISTORY:
249	! 07Jun02 - J.W. Larson <larson@mcs.anl.gov> - initial version
250	! ___________________________________________________________________
251
252	character(len=*),parameter :: myname_=myname//'::SpatialIntegralRAttrVDP_'
253
254	integer :: ierr, length
255	logical :: mySumWeights
256
257	! Argument Validity Checks
258
259	if(AttrVect_lsize(inAv) /= size(Weights)) then
260	ierr = AttrVect_lsize(inAv) - size(Weights)
261	write(stderr,'(3a,i8,a,i8)') myname_, &
262	':: inAv / Weights array length mismatch: ', &
263	' AttrVect_lsize(inAv) = ',AttrVect_lsize(inAv), &
264	' size(Weights) = ',size(Weights)
265	call die(myname_)
266	endif
267
268	if(present(SumWeights)) then
269	mySumWeights = SumWeights
270	if(.not. present(WeightTag)) then
271	write(stderr,'(3a)') myname_,':: FATAL--If the input argument SumWeights=.TRUE.,', &
272	' then the argument WeightTag must be provided.'
273	call die(myname_)
274	endif
275	else
276	mySumWeights = .FALSE.
277	endif
278
279	! Compute the sum
280
281	if(present(comm)) then ! compute distributed AllReduce-style sum:
282
283	if(mySumWeights) then ! return the spatial sum of the weights in outAV
284	call AttrVect_GlobalWeightedSumRAttr(inAV, outAV, Weights, &
285	comm, WeightTag)
286	else
287	call AttrVect_GlobalWeightedSumRAttr(inAV, outAV, Weights, comm)
288	endif
289
290	else ! compute local sum:
291
292	if(mySumWeights) then ! return the spatial sum of the weights in outAV
293	call AttrVect_LocalWeightedSumRAttr(inAV, outAV, Weights, &
294	WeightTag)
295	else
296	call AttrVect_LocalWeightedSumRAttr(inAV, outAV, Weights)
297	endif
298
299	endif ! if(present(comm))...
300
301	end subroutine SpatialIntegralRAttrVDP_
302
303
304	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
305	! Math and Computer Science Division, Argonne National Laboratory !
306	!BOP -------------------------------------------------------------------
307	!
308	! !IROUTINE: SpatialAverageRAttrVSP_ - Compute spatial average.
309	!
310	! !DESCRIPTION:
311	! This routine computes spatial averages of the {\tt REAL} attributes
312	! of the input {\tt AttrVect} argument {\tt inAv}.
313	! {\tt SpatialAverageRAttrV\_()} takes the input {\tt AttrVect} argument
314	! {\tt inAv} and computes the spatial average using weights
315	! stored in the {\tt REAL} array {\tt Weights}. The average of each
316	! {\tt REAL} attribute is returned in the output {\tt AttrVect} argument
317	! {\tt outAv}. If {\tt SpatialAverageRAttrV\_()} is invoked with the
318	! optional {\tt INTEGER} argument {\tt comm} (a Fortran MPI communicator
319	! handle), the summation operations for the average are completed on the
320	! local process, then reduced across the communicator, with all processes
321	! receiving the result.
322	!
323	! {\bf N.B.: } The local lengths of the {\tt AttrVect} argument {\tt inAv}
324	! and the input array {\tt Weights} must be equal. That is, there must
325	! be a one-to-one correspondence between the field point values stored
326	! in {\tt inAv} and the point weights stored in {\tt Weights}.
327	!
328	! {\bf N.B.: } The output {\tt AttrVect} argument {\tt outAv} is an
329	! allocated data structure. The user must deallocate it using the routine
330	! {\tt AttrVect\_clean()} when it is no longer needed. Failure to do so
331	! will result in a memory leak.
332	!
333	! !INTERFACE:
334
335	subroutine SpatialAverageRAttrVSP_(inAv, outAv, Weights, comm)
336
337	! ! USES:
338
339	use m_stdio
340	use m_die
341	use m_mpif90
342	use m_realkinds, only : SP, FP
343
344	use m_AttrVect, only : AttrVect
345	use m_AttrVect, only : AttrVect_init => init
346	use m_AttrVect, only : AttrVect_zero => zero
347	use m_AttrVect, only : AttrVect_clean => clean
348	use m_AttrVect, only : AttrVect_nRAttr => nRAttr
349	use m_AttrVect, only : AttrVect_indexRA => indexRA
350
351	use m_List, only : List
352	use m_List, only : List_nullify => nullify
353
354	implicit none
355
356	! !INPUT PARAMETERS:
357
358	type(AttrVect), intent(IN) :: inAv
359	real(SP), dimension(:), pointer :: Weights
360	integer, optional, intent(IN) :: comm
361
362	! !OUTPUT PARAMETERS:
363
364	type(AttrVect), intent(OUT) :: outAv
365
366	! !REVISION HISTORY:
367	! 10Jun02 - J.W. Larson <larson@mcs.anl.gov> - initial version
368	!EOP ___________________________________________________________________
369
370	character(len=*),parameter :: myname_=myname//'::SpatialAverageRAtttrVSP_'
371
372	type(AttrVect) :: integratedAv
373	type(List) :: nullIList
374	integer :: i, ierr, iweight
375
376	! Compute the spatial integral:
377
378	if(present(comm)) then
379	call SpatialIntegralV(inAv, integratedAv, Weights, &
380	.TRUE., 'weights', comm)
381	else
382	call SpatialIntegralV(inAv, integratedAv, Weights, .TRUE., 'weights')
383	endif
384
385	! Check value of summed weights (to avoid division by zero):
386
387	iweight = AttrVect_indexRA(integratedAv, 'weights')
388	if(integratedAv%rAttr(iweight, 1) == 0._FP) then
389	write(stderr,'(2a)') myname_, &
390	'::ERROR--Global sum of grid weights is zero.'
391	call die(myname_)
392	endif
393
394	! Initialize output AttrVect outAv:
395
396	call List_nullify(nullIList)
397	call AttrVect_init(outAv, iList=nullIList, rList=inAv%rList, lsize=1)
398	call AttrVect_zero(outAv)
399
400	! Divide by global weight sum to compute spatial averages from
401	! spatial integrals.
402
403	do i=1,AttrVect_nRAttr(outAv)
404	outAv%rAttr(i,1) = integratedAv%rAttr(i,1) &
405	/ integratedAv%rAttr(iweight,1)
406	end do
407
408	! Clean up temporary AttrVect:
409
410	call AttrVect_clean(integratedAv)
411
412	end subroutine SpatialAverageRAttrVSP_
413
414	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
415	! Math and Computer Science Division, Argonne National Laboratory !
416	! -------------------------------------------------------------------
417	!
418	! !IROUTINE: SpatialAverageRAttrVDP_ - Compute spatial average.
419	!
420	! !DESCRIPTION:
421	! Double pecision version of SpatialAverageRAttrVSP
422	!
423	! !INTERFACE:
424
425	subroutine SpatialAverageRAttrVDP_(inAv, outAv, Weights, comm)
426
427	! ! USES:
428
429	use m_stdio
430	use m_die
431	use m_mpif90
432	use m_realkinds, only : DP, FP
433
434	use m_AttrVect, only : AttrVect
435	use m_AttrVect, only : AttrVect_init => init
436	use m_AttrVect, only : AttrVect_zero => zero
437	use m_AttrVect, only : AttrVect_clean => clean
438	use m_AttrVect, only : AttrVect_nRAttr => nRAttr
439	use m_AttrVect, only : AttrVect_indexRA => indexRA
440
441	use m_List, only : List
442	use m_List, only : List_nullify => nullify
443
444	implicit none
445
446	! !INPUT PARAMETERS:
447
448	type(AttrVect), intent(IN) :: inAv
449	real(DP), dimension(:), pointer :: Weights
450	integer, optional, intent(IN) :: comm
451
452	! !OUTPUT PARAMETERS:
453
454	type(AttrVect), intent(OUT) :: outAv
455
456	! !REVISION HISTORY:
457	! 10Jun02 - J.W. Larson <larson@mcs.anl.gov> - initial version
458	! ___________________________________________________________________
459
460	character(len=*),parameter :: myname_=myname//'::SpatialAverageRAtttrVDP_'
461
462	type(AttrVect) :: integratedAv
463	type(List) :: nullIList
464	integer :: i, ierr, iweight
465
466	! Compute the spatial integral:
467
468	if(present(comm)) then
469	call SpatialIntegralV(inAv, integratedAv, Weights, &
470	.TRUE., 'weights', comm)
471	else
472	call SpatialIntegralV(inAv, integratedAv, Weights, .TRUE., 'weights')
473	endif
474
475	! Check value of summed weights (to avoid division by zero):
476
477	iweight = AttrVect_indexRA(integratedAv, 'weights')
478	if(integratedAv%rAttr(iweight, 1) == 0._FP) then
479	write(stderr,'(2a)') myname_, &
480	'::ERROR--Global sum of grid weights is zero.'
481	call die(myname_)
482	endif
483
484	! Initialize output AttrVect outAv:
485
486	call List_nullify(nullIList)
487	call AttrVect_init(outAv, iList=nullIList, rList=inAv%rList, lsize=1)
488	call AttrVect_zero(outAv)
489
490	! Divide by global weight sum to compute spatial averages from
491	! spatial integrals.
492
493	do i=1,AttrVect_nRAttr(outAv)
494	outAv%rAttr(i,1) = integratedAv%rAttr(i,1) &
495	/ integratedAv%rAttr(iweight,1)
496	end do
497
498	! Clean up temporary AttrVect:
499
500	call AttrVect_clean(integratedAv)
501
502	end subroutine SpatialAverageRAttrVDP_
503
504	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
505	! Math and Computer Science Division, Argonne National Laboratory !
506	!BOP -------------------------------------------------------------------
507	!
508	! !IROUTINE: MaskedSpatialIntegralRAttrVSP_ - Masked spatial integral.
509	!
510	! !DESCRIPTION:
511	! This routine computes masked spatial integrals of the {\tt REAL}
512	! attributes of the input {\tt AttrVect} argument {\tt inAv}, returning
513	! the masked integrals in the output {\tt AttrVect} argument {\tt outAv}.
514	! The masked integral is computed using weights stored in the input
515	! {\tt REAL} array argument {\tt SpatialWeights}. Integer masking (if
516	! desired) is provided in the optional input {\tt INTEGER} array {\tt iMask},
517	! and real masking (if desired) is provided in the optional input {\tt REAL}
518	! array {\tt rMask}. If {\tt SpatialIntegralRAttrV\_()} is invoked with the
519	! optional {\tt LOGICAL} input argument {\tt SumWeights} set as {\tt .TRUE.},
520	! then the weights are also summed and stored in {\tt outAv} (and can be
521	! referenced with the attribute name defined by the optional input
522	! {\tt CHARACTER} argument {\tt WeightSumTag}. If
523	! {\tt SpatialIntegralRAttrV\_()} is invoked with the optional {\tt INTEGER}
524	! argument {\tt comm} (a Fortran MPI communicator handle), the summation
525	! operations for the integral are completed on the local process, then
526	! reduced across the communicator, with all processes receiving the result.
527	! Otherwise, the integral is assumed to be local (or equivalent to a global
528	! address space).
529	!
530	! {\bf N.B.: } The local lengths of the {\tt AttrVect} argument {\tt inAv}
531	! and the input array {\tt Weights} must be equal. That is, there must be
532	! a one-to-one correspondence between the field point values stored
533	! in {\tt inAv} and the point weights stored in {\tt SpatialWeights}.
534	!
535	! {\bf N.B.: } If {\tt SpatialIntegralRAttrV\_()} is invoked with the
536	! optional {\tt LOGICAL} input argument {\tt SumWeights} set as {\tt .TRUE.}.
537	! In this case, the none of {\tt REAL} attribute tags in {\tt inAv} may be
538	! named the same as the string contained in {\tt WeightSumTag}, which is an
539	! attribute name reserved for the sum of the weights in the output {\tt AttrVect}
540	! {\tt outAv}.
541	!
542	! {\bf N.B.: } The output {\tt AttrVect} argument {\tt outAv} is an
543	! allocated data structure. The user must deallocate it using the routine
544	! {\tt AttrVect\_clean()} when it is no longer needed. Failure to do so
545	! will result in a memory leak.
546	!
547	! !INTERFACE:
548
549	subroutine MaskedSpatialIntegralRAttrVSP_(inAv, outAv, SpatialWeights, iMask, &
550	rMask, UseFastMethod, SumWeights, &
551	WeightSumTag, comm)
552
553	! ! USES:
554
555	use m_stdio
556	use m_die
557	use m_mpif90
558	use m_realkinds, only : SP, FP
559
560	use m_AttrVect, only : AttrVect
561	use m_AttrVect, only : AttrVect_lsize => lsize
562
563	use m_AttrVectReduce, only : AttrVect_GlobalWeightedSumRAttr => &
564	GlobalWeightedSumRAttr
565	use m_AttrVectReduce, only : AttrVect_LocalWeightedSumRAttr => &
566	LocalWeightedSumRAttr
567	implicit none
568
569	! !INPUT PARAMETERS:
570
571	type(AttrVect), intent(IN) :: inAv
572	real(SP),dimension(:), pointer :: SpatialWeights
573	integer, dimension(:), optional, pointer :: iMask
574	real(SP),dimension(:), optional, pointer :: rMask
575	logical, intent(IN) :: UseFastMethod
576	logical, optional, intent(IN) :: SumWeights
577	character(len=*), optional, intent(IN) :: WeightSumTag
578	integer, optional, intent(IN) :: comm
579
580	! !OUTPUT PARAMETERS:
581
582	type(AttrVect), intent(OUT) :: outAv
583
584	! !REVISION HISTORY:
585	! 10Jun02 - J.W. Larson <larson@mcs.anl.gov> - initial version
586	!
587	!EOP ___________________________________________________________________
588
589	character(len=*),parameter :: myname_=myname//'::MaskedSpatialIntegralRAttrVSP_'
590
591	integer :: i, ierr, length
592	logical :: mySumWeights
593	real(FP), dimension(:), pointer :: Weights
594
595	! Argument Validity Checks
596
597	if(AttrVect_lsize(inAv) /= size(SpatialWeights)) then
598	ierr = AttrVect_lsize(inAv) - size(SpatialWeights)
599	write(stderr,'(3a,i8,a,i8)') myname_, &
600	':: inAv / SpatialWeights array length mismatch: ', &
601	' AttrVect_lsize(inAv) = ',AttrVect_lsize(inAv), &
602	' size(SpatialWeights) = ',size(SpatialWeights)
603	call die(myname_)
604	endif
605
606	if(present(iMask)) then ! make sure it is the right length
607	if(AttrVect_lsize(inAv) /= size(iMask)) then
608	ierr = AttrVect_lsize(inAv) - size(iMask)
609	write(stderr,'(3a,i8,a,i8)') myname_, &
610	':: inAv / iMask array length mismatch: ', &
611	' AttrVect_lsize(inAv) = ',AttrVect_lsize(inAv), &
612	' size(iMask) = ',size(iMask)
613	call die(myname_)
614	endif
615	endif
616
617	if(present(rMask)) then ! make sure it is the right length
618	if(AttrVect_lsize(inAv) /= size(rMask)) then
619	ierr = AttrVect_lsize(inAv) - size(rMask)
620	write(stderr,'(3a,i8,a,i8)') myname_, &
621	':: inAv / rMask array length mismatch: ', &
622	' AttrVect_lsize(inAv) = ',AttrVect_lsize(inAv), &
623	' size(rMask) = ',size(rMask)
624	call die(myname_)
625	endif
626	endif
627
628	if(present(SumWeights)) then
629	mySumWeights = SumWeights
630	if(.not. present(WeightSumTag)) then
631	write(stderr,'(3a)') myname_,':: FATAL--If the input argument SumWeights=.TRUE.,', &
632	' then the argument WeightSumTag must be provided.'
633	call die(myname_)
634	endif
635	else
636	mySumWeights = .FALSE.
637	endif
638
639	! Create a common Weights(:) array...
640
641	length = AttrVect_lsize(inAv)
642
643	allocate(Weights(length), stat=ierr)
644	if(ierr /= 0) then
645	write(stderr,'(3a,i8)') myname_,':: allocate(Weights(...) failed,', &
646	' ierr=',ierr
647	call die(myname_)
648	endif
649
650	! Combine weights and masks into a common Weights(:) array...
651
652	if(UseFastMethod) then ! form the product of iMask, rMask, and SpatialWeights
653
654	if(present(rMask)) then ! use it to form Weights(:)
655	if(present(iMask)) then ! use it and rMask to form Weights(:)
656	do i=1,length
657	Weights(i) = rMask(i) * SpatialWeights(i) * iMask(i)
658	end do
659	else
660	do i=1,length
661	Weights(i) = rMask(i) * SpatialWeights(i)
662	end do
663	endif ! if(present(iMask))...
664	else
665	if(present(iMask)) then
666	do i=1,length
667	Weights(i) = SpatialWeights(i) * iMask(i)
668	end do
669	else
670	do i=1,length
671	Weights(i) = SpatialWeights(i)
672	end do
673	endif ! if(present(iMask))...
674	endif ! if(present(rMask))...
675
676
677	else ! Scan iMask and rMask carefully and set Weights(i) to zero
678	! when iMask(i) or rMask(i) is zero. This avoids round-off
679	! effects from products and promotion of integers to reals.
680
681	if(present(rMask)) then ! use it to form Weights(:)
682	if(present(iMask)) then ! use it and rMask to form Weights(:)
683	do i=1,length
684	select case(iMask(i))
685	case(0)
686	Weights(i) = 0._FP
687	case(1)
688	if(rMask(i) == 1._FP) then
689	Weights(i) = SpatialWeights(i)
690	elseif(rMask(i) == 0._FP) then
691	Weights(i) = 0._FP
692	elseif((rMask(i) > 0._FP) .and. (rMask(i) < 1._FP)) then
693	Weights(i) = rMask(i) * SpatialWeights(i)
694	else ! rMask(i) < 0. or rMask(i) > 1.
695	write(stderr,'(3a,i8,a,f10.7)') myname_, &
696	':: invalid value for real', &
697	'mask entry rMask(',i,') = ',rMask(i)
698	call die(myname_)
699	endif
700	case default
701	write(stderr,'(3a,i8,a,i8)') myname_, &
702	':: invalid value for integer', &
703	'mask entry iMask(',i,') = ',iMask(i)
704	call die(myname_)
705	end select
706	end do
707	else
708	do i=1,length
709	if(rMask(i) == 1._FP) then
710	Weights(i) = SpatialWeights(i)
711	elseif(rMask(i) == 0._FP) then
712	Weights(i) = 0._FP
713	elseif((rMask(i) > 0._FP) .and. (rMask(i) < 1._FP)) then
714	Weights(i) = rMask(i) * SpatialWeights(i)
715	else ! rMask(i) < 0. or rMask(i) > 1.
716	write(stderr,'(3a,i8,a,e10.6)') myname_, &
717	':: invalid value for real', &
718	'mask entry rMask(',i,') = ',rMask(i)
719	call die(myname_)
720	endif
721	end do
722	endif ! if(present(iMask))...
723	else ! no rMask present...
724	if(present(iMask)) then ! check iMask entries...
725	do i=1,length
726	select case(iMask(i))
727	case(0)
728	Weights(i) = 0._FP
729	case(1)
730	Weights(i) = SpatialWeights(i)
731	case default
732	write(stderr,'(3a,i8,a,i8)') myname_, &
733	':: invalid value for integer', &
734	'mask entry iMask(',i,') = ',iMask(i)
735	call die(myname_)
736	end select
737	end do
738	else ! straight assignment of SpatialWeights(:)
739	do i=1,length
740	Weights(i) = SpatialWeights(i)
741	end do
742	endif ! if(present(iMask))...
743	endif ! if(present(rMask))...
744
745
746	endif ! if(UseFastMethod)
747
748	! Now that the weights are combined into a common Weights(:),
749	! compute the masked weighted sum:
750
751	if(present(comm)) then ! compute distributed AllReduce-style sum:
752
753	if(mySumWeights) then ! return the global sum of the weights in outAV
754	call AttrVect_GlobalWeightedSumRAttr(inAV, outAV, Weights, &
755	comm, WeightSumTag)
756	else
757	call AttrVect_GlobalWeightedSumRAttr(inAV, outAV, Weights, comm)
758	endif
759
760	else ! compute local sum:
761
762	if(mySumWeights) then ! return the global sum of the weights in outAV
763	call AttrVect_LocalWeightedSumRAttr(inAV, outAV, Weights, &
764	WeightSumAttr=WeightSumTag)
765	else
766	call AttrVect_LocalWeightedSumRAttr(inAV, outAV, Weights)
767	endif
768
769	endif ! if(present(comm))...
770
771	! Clean up the allocated Weights(:) array
772
773	deallocate(Weights, stat=ierr)
774	if(ierr /= 0) then
775	write(stderr,'(3a,i8)') myname_,':: deallocate(Weights(...) failed,', &
776	' ierr=',ierr
777	call die(myname_)
778	endif
779
780	end subroutine MaskedSpatialIntegralRAttrVSP_
781
782	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
783	! Math and Computer Science Division, Argonne National Laboratory !
784	! -------------------------------------------------------------------
785	!
786	! !IROUTINE: MaskedSpatialIntegralRAttrVDP_ - Masked spatial integral.
787	!
788	! !DESCRIPTION:
789	! Double precision version of MaskedSpatialIntegralRAttrVSP_
790	!
791	! !INTERFACE:
792
793	subroutine MaskedSpatialIntegralRAttrVDP_(inAv, outAv, SpatialWeights, iMask, &
794	rMask, UseFastMethod, SumWeights, &
795	WeightSumTag, comm)
796
797	! ! USES:
798
799	use m_stdio
800	use m_die
801	use m_mpif90
802	use m_realkinds, only : DP, FP
803
804	use m_AttrVect, only : AttrVect
805	use m_AttrVect, only : AttrVect_lsize => lsize
806
807	use m_AttrVectReduce, only : AttrVect_GlobalWeightedSumRAttr => &
808	GlobalWeightedSumRAttr
809	use m_AttrVectReduce, only : AttrVect_LocalWeightedSumRAttr => &
810	LocalWeightedSumRAttr
811	implicit none
812
813	! !INPUT PARAMETERS:
814
815	type(AttrVect), intent(IN) :: inAv
816	real(DP),dimension(:), pointer :: SpatialWeights
817	integer, dimension(:), optional, pointer :: iMask
818	real(DP),dimension(:), optional, pointer :: rMask
819	logical, intent(IN) :: UseFastMethod
820	logical, optional, intent(IN) :: SumWeights
821	character(len=*), optional, intent(IN) :: WeightSumTag
822	integer, optional, intent(IN) :: comm
823
824	! !OUTPUT PARAMETERS:
825
826	type(AttrVect), intent(OUT) :: outAv
827
828	! !REVISION HISTORY:
829	! 10Jun02 - J.W. Larson <larson@mcs.anl.gov> - initial version
830	! ___________________________________________________________________
831
832	character(len=*),parameter :: myname_=myname//'::MaskedSpatialIntegralRAttrVDP_'
833
834	integer :: i, ierr, length
835	logical :: mySumWeights
836	real(FP), dimension(:), pointer :: Weights
837
838	! Argument Validity Checks
839
840	if(AttrVect_lsize(inAv) /= size(SpatialWeights)) then
841	ierr = AttrVect_lsize(inAv) - size(SpatialWeights)
842	write(stderr,'(3a,i8,a,i8)') myname_, &
843	':: inAv / SpatialWeights array length mismatch: ', &
844	' AttrVect_lsize(inAv) = ',AttrVect_lsize(inAv), &
845	' size(SpatialWeights) = ',size(SpatialWeights)
846	call die(myname_)
847	endif
848
849	if(present(iMask)) then ! make sure it is the right length
850	if(AttrVect_lsize(inAv) /= size(iMask)) then
851	ierr = AttrVect_lsize(inAv) - size(iMask)
852	write(stderr,'(3a,i8,a,i8)') myname_, &
853	':: inAv / iMask array length mismatch: ', &
854	' AttrVect_lsize(inAv) = ',AttrVect_lsize(inAv), &
855	' size(iMask) = ',size(iMask)
856	call die(myname_)
857	endif
858	endif
859
860	if(present(rMask)) then ! make sure it is the right length
861	if(AttrVect_lsize(inAv) /= size(rMask)) then
862	ierr = AttrVect_lsize(inAv) - size(rMask)
863	write(stderr,'(3a,i8,a,i8)') myname_, &
864	':: inAv / rMask array length mismatch: ', &
865	' AttrVect_lsize(inAv) = ',AttrVect_lsize(inAv), &
866	' size(rMask) = ',size(rMask)
867	call die(myname_)
868	endif
869	endif
870
871	if(present(SumWeights)) then
872	mySumWeights = SumWeights
873	if(.not. present(WeightSumTag)) then
874	write(stderr,'(3a)') myname_,':: FATAL--If the input argument SumWeights=.TRUE.,', &
875	' then the argument WeightSumTag must be provided.'
876	call die(myname_)
877	endif
878	else
879	mySumWeights = .FALSE.
880	endif
881
882	! Create a common Weights(:) array...
883
884	length = AttrVect_lsize(inAv)
885
886	allocate(Weights(length), stat=ierr)
887	if(ierr /= 0) then
888	write(stderr,'(3a,i8)') myname_,':: allocate(Weights(...) failed,', &
889	' ierr=',ierr
890	call die(myname_)
891	endif
892
893	! Combine weights and masks into a common Weights(:) array...
894
895	if(UseFastMethod) then ! form the product of iMask, rMask, and SpatialWeights
896
897	if(present(rMask)) then ! use it to form Weights(:)
898	if(present(iMask)) then ! use it and rMask to form Weights(:)
899	do i=1,length
900	Weights(i) = rMask(i) * SpatialWeights(i) * iMask(i)
901	end do
902	else
903	do i=1,length
904	Weights(i) = rMask(i) * SpatialWeights(i)
905	end do
906	endif ! if(present(iMask))...
907	else
908	if(present(iMask)) then
909	do i=1,length
910	Weights(i) = SpatialWeights(i) * iMask(i)
911	end do
912	else
913	do i=1,length
914	Weights(i) = SpatialWeights(i)
915	end do
916	endif ! if(present(iMask))...
917	endif ! if(present(rMask))...
918
919
920	else ! Scan iMask and rMask carefully and set Weights(i) to zero
921	! when iMask(i) or rMask(i) is zero. This avoids round-off
922	! effects from products and promotion of integers to reals.
923
924	if(present(rMask)) then ! use it to form Weights(:)
925	if(present(iMask)) then ! use it and rMask to form Weights(:)
926	do i=1,length
927	select case(iMask(i))
928	case(0)
929	Weights(i) = 0._FP
930	case(1)
931	if(rMask(i) == 1._FP) then
932	Weights(i) = SpatialWeights(i)
933	elseif(rMask(i) == 0._FP) then
934	Weights(i) = 0._FP
935	elseif((rMask(i) > 0._FP) .and. (rMask(i) < 1._FP)) then
936	Weights(i) = rMask(i) * SpatialWeights(i)
937	else ! rMask(i) < 0. or rMask(i) > 1.
938	write(stderr,'(3a,i8,a,f10.7)') myname_, &
939	':: invalid value for real', &
940	'mask entry rMask(',i,') = ',rMask(i)
941	call die(myname_)
942	endif
943	case default
944	write(stderr,'(3a,i8,a,i8)') myname_, &
945	':: invalid value for integer', &
946	'mask entry iMask(',i,') = ',iMask(i)
947	call die(myname_)
948	end select
949	end do
950	else
951	do i=1,length
952	if(rMask(i) == 1._FP) then
953	Weights(i) = SpatialWeights(i)
954	elseif(rMask(i) == 0._FP) then
955	Weights(i) = 0._FP
956	elseif((rMask(i) > 0._FP) .and. (rMask(i) < 1._FP)) then
957	Weights(i) = rMask(i) * SpatialWeights(i)
958	else ! rMask(i) < 0. or rMask(i) > 1.
959	write(stderr,'(3a,i8,a,e10.6)') myname_, &
960	':: invalid value for real', &
961	'mask entry rMask(',i,') = ',rMask(i)
962	call die(myname_)
963	endif
964	end do
965	endif ! if(present(iMask))...
966	else ! no rMask present...
967	if(present(iMask)) then ! check iMask entries...
968	do i=1,length
969	select case(iMask(i))
970	case(0)
971	Weights(i) = 0._FP
972	case(1)
973	Weights(i) = SpatialWeights(i)
974	case default
975	write(stderr,'(3a,i8,a,i8)') myname_, &
976	':: invalid value for integer', &
977	'mask entry iMask(',i,') = ',iMask(i)
978	call die(myname_)
979	end select
980	end do
981	else ! straight assignment of SpatialWeights(:)
982	do i=1,length
983	Weights(i) = SpatialWeights(i)
984	end do
985	endif ! if(present(iMask))...
986	endif ! if(present(rMask))...
987
988
989	endif ! if(UseFastMethod)
990
991	! Now that the weights are combined into a common Weights(:),
992	! compute the masked weighted sum:
993
994	if(present(comm)) then ! compute distributed AllReduce-style sum:
995
996	if(mySumWeights) then ! return the global sum of the weights in outAV
997	call AttrVect_GlobalWeightedSumRAttr(inAV, outAV, Weights, &
998	comm, WeightSumTag)
999	else
1000	call AttrVect_GlobalWeightedSumRAttr(inAV, outAV, Weights, comm)
1001	endif
1002
1003	else ! compute local sum:
1004
1005	if(mySumWeights) then ! return the global sum of the weights in outAV
1006	call AttrVect_LocalWeightedSumRAttr(inAV, outAV, Weights, &
1007	WeightSumAttr=WeightSumTag)
1008	else
1009	call AttrVect_LocalWeightedSumRAttr(inAV, outAV, Weights)
1010	endif
1011
1012	endif ! if(present(comm))...
1013
1014	! Clean up the allocated Weights(:) array
1015
1016	deallocate(Weights, stat=ierr)
1017	if(ierr /= 0) then
1018	write(stderr,'(3a,i8)') myname_,':: deallocate(Weights(...) failed,', &
1019	' ierr=',ierr
1020	call die(myname_)
1021	endif
1022
1023	end subroutine MaskedSpatialIntegralRAttrVDP_
1024
1025
1026	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1027	! Math and Computer Science Division, Argonne National Laboratory !
1028	!BOP -------------------------------------------------------------------
1029	!
1030	! !IROUTINE: MaskedSpatialAverageRAttrVSP_ - Masked spatial average.
1031	!
1032	! !DESCRIPTION: [NEEDS LOTS of work...]
1033	! This routine computes spatial integrals of the {\tt REAL} attributes
1034	! of the {\tt REAL} attributes of the input {\tt AttrVect} argument
1035	! {\tt inAv}. {\tt SpatialIntegralRAttrV\_()} takes the input
1036	! {\tt AttrVect} argument {\tt inAv} and computes the spatial
1037	! integral using weights stored in the input {\tt REAL} array argument
1038	! {\tt Weights}. The integral of each {\tt REAL} attribute is returned
1039	! in the output {\tt AttrVect} argument {\tt outAv}. If
1040	! {\tt SpatialIntegralRAttrV\_()} is invoked with the optional {\tt LOGICAL}
1041	! input argument {\tt SumWeights} set as {\tt .TRUE.}, then the weights
1042	! are also summed and stored in {\tt outAv} (and can be referenced with
1043	! the attribute name {\tt WeightTag}. If {\tt SpatialIntegralRAttrV\_()} is
1044	! invoked with the optional {\tt INTEGER} argument {\tt comm} (a Fortran
1045	! MPI communicator handle), the summation operations for the integral are
1046	! completed on the local process, then reduced across the communicator,
1047	! with all processes receiving the result.
1048	!
1049	! {\bf N.B.: } The local lengths of the {\tt AttrVect} argument {\tt inAv}
1050	! and the input array {\tt Weights} must be equal. That is, there must be
1051	! a one-to-one correspondence between the field point values stored
1052	! in {\tt inAv} and the point weights stored in {\tt Weights}.
1053	!
1054	! {\bf N.B.: } If {\tt SpatialIntegralRAttrV\_()} is invoked with the
1055	! optional {\tt LOGICAL} input argument {\tt SumWeights} set as {\tt .TRUE.}.
1056	! In this case, the none of {\tt REAL} attribute tags in {\tt inAv} may be
1057	! named the same as the string contained in {\tt WeightTag}, which is an
1058	! attribute name reserved for the sum of the weights in the output {\tt AttrVect}
1059	! {\tt outAv}.
1060	!
1061	! {\bf N.B.: } The output {\tt AttrVect} argument {\tt outAv} is an
1062	! allocated data structure. The user must deallocate it using the routine
1063	! {\tt AttrVect\_clean()} when it is no longer needed. Failure to do so
1064	! will result in a memory leak.
1065	!
1066	! !INTERFACE:
1067
1068	subroutine MaskedSpatialAverageRAttrVSP_(inAv, outAv, SpatialWeights, iMask, &
1069	rMask, UseFastMethod, comm)
1070
1071	! ! USES:
1072
1073	use m_stdio
1074	use m_die
1075	use m_mpif90
1076	use m_realkinds, only : SP, FP
1077
1078	use m_AttrVect, only : AttrVect
1079	use m_AttrVect, only : AttrVect_init => init
1080	use m_AttrVect, only : AttrVect_zero => zero
1081	use m_AttrVect, only : AttrVect_clean => clean
1082	use m_AttrVect, only : AttrVect_lsize => lsize
1083	use m_AttrVect, only : AttrVect_nRAttr => nRAttr
1084	use m_AttrVect, only : AttrVect_indexRA => indexRA
1085
1086	use m_List, only : List
1087	use m_List, only : List_nullify => nullify
1088
1089	implicit none
1090
1091	! !INPUT PARAMETERS:
1092
1093	type(AttrVect), intent(IN) :: inAv
1094	real(SP), dimension(:), pointer :: SpatialWeights
1095	integer, dimension(:), optional, pointer :: iMask
1096	real(SP),dimension(:), optional, pointer :: rMask
1097	logical, intent(IN) :: UseFastMethod
1098	integer, optional, intent(IN) :: comm
1099
1100	! !OUTPUT PARAMETERS:
1101
1102	type(AttrVect), intent(OUT) :: outAv
1103
1104	! !REVISION HISTORY:
1105	! 11Jun02 - J.W. Larson <larson@mcs.anl.gov> - initial version
1106	!EOP ___________________________________________________________________
1107
1108	character(len=*),parameter :: myname_=myname//'::MaskedSpatialAverageRAttrVSP_'
1109
1110	type(AttrVect) :: integratedAv
1111	type(List) :: nullIList
1112
1113	integer :: i, ierr, length, iweight
1114	logical :: mySumWeights
1115
1116	! Argument Validity Checks
1117
1118	if(AttrVect_lsize(inAv) /= size(SpatialWeights)) then
1119	ierr = AttrVect_lsize(inAv) - size(SpatialWeights)
1120	write(stderr,'(3a,i8,a,i8)') myname_, &
1121	':: inAv / SpatialWeights array length mismatch: ', &
1122	' AttrVect_lsize(inAv) = ',AttrVect_lsize(inAv), &
1123	' size(SpatialWeights) = ',size(SpatialWeights)
1124	call die(myname_)
1125	endif
1126
1127	if(present(iMask)) then ! make sure it is the right length
1128	if(AttrVect_lsize(inAv) /= size(iMask)) then
1129	ierr = AttrVect_lsize(inAv) - size(iMask)
1130	write(stderr,'(3a,i8,a,i8)') myname_, &
1131	':: inAv / iMask array length mismatch: ', &
1132	' AttrVect_lsize(inAv) = ',AttrVect_lsize(inAv), &
1133	' size(iMask) = ',size(iMask)
1134	call die(myname_)
1135	endif
1136	endif
1137
1138	if(present(rMask)) then ! make sure it is the right length
1139	if(AttrVect_lsize(inAv) /= size(rMask)) then
1140	ierr = AttrVect_lsize(inAv) - size(rMask)
1141	write(stderr,'(3a,i8,a,i8)') myname_, &
1142	':: inAv / rMask array length mismatch: ', &
1143	' AttrVect_lsize(inAv) = ',AttrVect_lsize(inAv), &
1144	' size(rMask) = ',size(rMask)
1145	call die(myname_)
1146	endif
1147	endif
1148
1149	! Compute the masked weighted sum, including the sum of the
1150	! masked weights.
1151
1152	if(present(comm)) then ! communicator handle present
1153
1154	if(present(iMask)) then
1155
1156	if(present(rMask)) then
1157	call MaskedSpatialIntegralV(inAv, integratedAv, SpatialWeights, &
1158	iMask, rMask, UseFastMethod, .TRUE., &
1159	'MaskedWeightsSum', comm)
1160	else ! no rMask
1161	call MaskedSpatialIntegralV(inAv, integratedAv, SpatialWeights, &
1162	iMask=iMask, UseFastMethod=UseFastMethod, &
1163	SumWeights=.TRUE., &
1164	WeightSumTag='MaskedWeightsSum', &
1165	comm=comm)
1166	endif ! if(present(rMask))...
1167
1168	else ! no iMask present...
1169
1170	if(present(rMask)) then
1171	call MaskedSpatialIntegralV(inAv, integratedAv, SpatialWeights, &
1172	rMask=rMask, UseFastMethod=UseFastMethod, &
1173	SumWeights=.TRUE., &
1174	WeightSumTag='MaskedWeightsSum', &
1175	comm=comm)
1176	else ! neither rMask nor iMask present:
1177	call MaskedSpatialIntegralV(inAv, integratedAv, SpatialWeights, &
1178	UseFastMethod=UseFastMethod, &
1179	SumWeights=.TRUE., &
1180	WeightSumTag='MaskedWeightsSum', &
1181	comm=comm)
1182	endif ! if(present(rMask))...
1183
1184	endif ! if(present(iMask))...
1185
1186	else ! no communicator handle present
1187
1188	if(present(iMask)) then
1189
1190	if(present(rMask)) then
1191	call MaskedSpatialIntegralV(inAv, integratedAv, SpatialWeights, &
1192	iMask, rMask, UseFastMethod, .TRUE., &
1193	'MaskedWeightsSum')
1194	else ! no rMask
1195	call MaskedSpatialIntegralV(inAv, integratedAv, SpatialWeights, &
1196	iMask=iMask, UseFastMethod=UseFastMethod, &
1197	SumWeights=.TRUE., &
1198	WeightSumTag='MaskedWeightsSum')
1199	endif ! if(present(rMask))...
1200
1201	else ! no iMask present...
1202
1203	if(present(rMask)) then
1204	call MaskedSpatialIntegralV(inAv, integratedAv, SpatialWeights, &
1205	rMask=rMask, UseFastMethod=UseFastMethod, &
1206	SumWeights=.TRUE., &
1207	WeightSumTag='MaskedWeightsSum')
1208	else ! neither rMask nor iMask present:
1209	call MaskedSpatialIntegralV(inAv, integratedAv, SpatialWeights, &
1210	UseFastMethod=UseFastMethod, &
1211	SumWeights=.TRUE., &
1212	WeightSumTag='MaskedWeightsSum')
1213	endif ! if(present(rMask))...
1214
1215	endif ! if(present(iMask))...
1216
1217	endif ! if(present(comm))...
1218
1219	! At this point, integratedAv containes the masked spatial integrals
1220	! of the REAL attributes of inAv, along with the sum of the weights.
1221	! to compute the masked spatial average
1222
1223	! Check value of summed weights (to avoid division by zero):
1224
1225	iweight = AttrVect_indexRA(integratedAv, 'MaskedWeightsSum')
1226	if(integratedAv%rAttr(iweight, 1) == 0._FP) then
1227	write(stderr,'(2a)') myname_, &
1228	'::ERROR--Global sum of grid weights is zero.'
1229	call die(myname_)
1230	endif
1231
1232	! Initialize output AttrVect outAv:
1233
1234	call List_nullify(nullIList)
1235	call AttrVect_init(outAv, iList=nullIList, rList=inAv%rList, lsize=1)
1236	call AttrVect_zero(outAv)
1237
1238	! Divide by global weight sum to compute spatial averages from
1239	! spatial integrals.
1240
1241	do i=1,AttrVect_nRAttr(outAv)
1242	outAv%rAttr(i,1) = integratedAv%rAttr(i,1) &
1243	/ integratedAv%rAttr(iweight,1)
1244	end do
1245
1246	! Clean up temporary AttrVect:
1247
1248	call AttrVect_clean(integratedAv)
1249
1250	end subroutine MaskedSpatialAverageRAttrVSP_
1251
1252	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1253	! Math and Computer Science Division, Argonne National Laboratory !
1254	! -------------------------------------------------------------------
1255	!
1256	! !IROUTINE: MaskedSpatialAverageRAttrVDP_ - Masked spatial average.
1257	!
1258	! !DESCRIPTION: [NEEDS LOTS of work...]
1259	! Double precision interface version of MaskedSpatialAverageRAttrVSP_.
1260	!
1261	! !INTERFACE:
1262
1263	subroutine MaskedSpatialAverageRAttrVDP_(inAv, outAv, SpatialWeights, iMask, &
1264	rMask, UseFastMethod, comm)
1265
1266	! ! USES:
1267
1268	use m_stdio
1269	use m_die
1270	use m_mpif90
1271	use m_realkinds, only : DP, FP
1272
1273	use m_AttrVect, only : AttrVect
1274	use m_AttrVect, only : AttrVect_init => init
1275	use m_AttrVect, only : AttrVect_zero => zero
1276	use m_AttrVect, only : AttrVect_clean => clean
1277	use m_AttrVect, only : AttrVect_lsize => lsize
1278	use m_AttrVect, only : AttrVect_nRAttr => nRAttr
1279	use m_AttrVect, only : AttrVect_indexRA => indexRA
1280
1281	use m_List, only : List
1282	use m_List, only : List_nullify => nullify
1283
1284	implicit none
1285
1286	! !INPUT PARAMETERS:
1287
1288	type(AttrVect), intent(IN) :: inAv
1289	real(DP), dimension(:), pointer :: SpatialWeights
1290	integer, dimension(:), optional, pointer :: iMask
1291	real(DP),dimension(:), optional, pointer :: rMask
1292	logical, intent(IN) :: UseFastMethod
1293	integer, optional, intent(IN) :: comm
1294
1295	! !OUTPUT PARAMETERS:
1296
1297	type(AttrVect), intent(OUT) :: outAv
1298
1299	! !REVISION HISTORY:
1300	! 11Jun02 - J.W. Larson <larson@mcs.anl.gov> - initial version
1301	! ___________________________________________________________________
1302
1303	character(len=*),parameter :: myname_=myname//'::MaskedSpatialAverageRAttrVDP_'
1304
1305	type(AttrVect) :: integratedAv
1306	type(List) :: nullIList
1307
1308	integer :: i, ierr, length, iweight
1309	logical :: mySumWeights
1310
1311	! Argument Validity Checks
1312
1313	if(AttrVect_lsize(inAv) /= size(SpatialWeights)) then
1314	ierr = AttrVect_lsize(inAv) - size(SpatialWeights)
1315	write(stderr,'(3a,i8,a,i8)') myname_, &
1316	':: inAv / SpatialWeights array length mismatch: ', &
1317	' AttrVect_lsize(inAv) = ',AttrVect_lsize(inAv), &
1318	' size(SpatialWeights) = ',size(SpatialWeights)
1319	call die(myname_)
1320	endif
1321
1322	if(present(iMask)) then ! make sure it is the right length
1323	if(AttrVect_lsize(inAv) /= size(iMask)) then
1324	ierr = AttrVect_lsize(inAv) - size(iMask)
1325	write(stderr,'(3a,i8,a,i8)') myname_, &
1326	':: inAv / iMask array length mismatch: ', &
1327	' AttrVect_lsize(inAv) = ',AttrVect_lsize(inAv), &
1328	' size(iMask) = ',size(iMask)
1329	call die(myname_)
1330	endif
1331	endif
1332
1333	if(present(rMask)) then ! make sure it is the right length
1334	if(AttrVect_lsize(inAv) /= size(rMask)) then
1335	ierr = AttrVect_lsize(inAv) - size(rMask)
1336	write(stderr,'(3a,i8,a,i8)') myname_, &
1337	':: inAv / rMask array length mismatch: ', &
1338	' AttrVect_lsize(inAv) = ',AttrVect_lsize(inAv), &
1339	' size(rMask) = ',size(rMask)
1340	call die(myname_)
1341	endif
1342	endif
1343
1344	! Compute the masked weighted sum, including the sum of the
1345	! masked weights.
1346
1347	if(present(comm)) then ! communicator handle present
1348
1349	if(present(iMask)) then
1350
1351	if(present(rMask)) then
1352	call MaskedSpatialIntegralV(inAv, integratedAv, SpatialWeights, &
1353	iMask, rMask, UseFastMethod, .TRUE., &
1354	'MaskedWeightsSum', comm)
1355	else ! no rMask
1356	call MaskedSpatialIntegralV(inAv, integratedAv, SpatialWeights, &
1357	iMask=iMask, UseFastMethod=UseFastMethod, &
1358	SumWeights=.TRUE., &
1359	WeightSumTag='MaskedWeightsSum', &
1360	comm=comm)
1361	endif ! if(present(rMask))...
1362
1363	else ! no iMask present...
1364
1365	if(present(rMask)) then
1366	call MaskedSpatialIntegralV(inAv, integratedAv, SpatialWeights, &
1367	rMask=rMask, UseFastMethod=UseFastMethod, &
1368	SumWeights=.TRUE., &
1369	WeightSumTag='MaskedWeightsSum', &
1370	comm=comm)
1371	else ! neither rMask nor iMask present:
1372	call MaskedSpatialIntegralV(inAv, integratedAv, SpatialWeights, &
1373	UseFastMethod=UseFastMethod, &
1374	SumWeights=.TRUE., &
1375	WeightSumTag='MaskedWeightsSum', &
1376	comm=comm)
1377	endif ! if(present(rMask))...
1378
1379	endif ! if(present(iMask))...
1380
1381	else ! no communicator handle present
1382
1383	if(present(iMask)) then
1384
1385	if(present(rMask)) then
1386	call MaskedSpatialIntegralV(inAv, integratedAv, SpatialWeights, &
1387	iMask, rMask, UseFastMethod, .TRUE., &
1388	'MaskedWeightsSum')
1389	else ! no rMask
1390	call MaskedSpatialIntegralV(inAv, integratedAv, SpatialWeights, &
1391	iMask=iMask, UseFastMethod=UseFastMethod, &
1392	SumWeights=.TRUE., &
1393	WeightSumTag='MaskedWeightsSum')
1394	endif ! if(present(rMask))...
1395
1396	else ! no iMask present...
1397
1398	if(present(rMask)) then
1399	call MaskedSpatialIntegralV(inAv, integratedAv, SpatialWeights, &
1400	rMask=rMask, UseFastMethod=UseFastMethod, &
1401	SumWeights=.TRUE., &
1402	WeightSumTag='MaskedWeightsSum')
1403	else ! neither rMask nor iMask present:
1404	call MaskedSpatialIntegralV(inAv, integratedAv, SpatialWeights, &
1405	UseFastMethod=UseFastMethod, &
1406	SumWeights=.TRUE., &
1407	WeightSumTag='MaskedWeightsSum')
1408	endif ! if(present(rMask))...
1409
1410	endif ! if(present(iMask))...
1411
1412	endif ! if(present(comm))...
1413
1414	! At this point, integratedAv containes the masked spatial integrals
1415	! of the REAL attributes of inAv, along with the sum of the weights.
1416	! to compute the masked spatial average
1417
1418	! Check value of summed weights (to avoid division by zero):
1419
1420	iweight = AttrVect_indexRA(integratedAv, 'MaskedWeightsSum')
1421	if(integratedAv%rAttr(iweight, 1) == 0._FP) then
1422	write(stderr,'(2a)') myname_, &
1423	'::ERROR--Global sum of grid weights is zero.'
1424	call die(myname_)
1425	endif
1426
1427	! Initialize output AttrVect outAv:
1428
1429	call List_nullify(nullIList)
1430	call AttrVect_init(outAv, iList=nullIList, rList=inAv%rList, lsize=1)
1431	call AttrVect_zero(outAv)
1432
1433	! Divide by global weight sum to compute spatial averages from
1434	! spatial integrals.
1435
1436	do i=1,AttrVect_nRAttr(outAv)
1437	outAv%rAttr(i,1) = integratedAv%rAttr(i,1) &
1438	/ integratedAv%rAttr(iweight,1)
1439	end do
1440
1441	! Clean up temporary AttrVect:
1442
1443	call AttrVect_clean(integratedAv)
1444
1445	end subroutine MaskedSpatialAverageRAttrVDP_
1446
1447	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1448	! Math and Computer Science Division, Argonne National Laboratory !
1449	!BOP -------------------------------------------------------------------
1450	!
1451	! !IROUTINE: PairedSpatialIntegralRAttrVSP_ - Do two spatial integrals at once.
1452	!
1453	! !DESCRIPTION:
1454	! This routine computes spatial integrals of the {\tt REAL} attributes
1455	! of the {\tt REAL} attributes of the input {\tt AttrVect} arguments
1456	! {\tt inAv1} and {\tt inAv2}, returning the integrals in the output
1457	! {\tt AttrVect} arguments {\tt outAv1} and {\tt outAv2}, respectively .
1458	! The integrals of {\tt inAv1} and {\tt inAv2} are computed using
1459	! spatial weights stored in the input {\tt REAL} array arguments
1460	! {\tt Weights1} and {\tt Weights2}, respectively.
1461	! If {\tt SpatialIntegralRAttrV\_()} is invoked with the optional
1462	! {\tt LOGICAL} input argument
1463	! {\tt SumWeights} set as {\tt .TRUE.}, then the weights are also summed
1464	! and stored in {\tt outAv1} and {\tt outAv2}, and can be referenced with
1465	! the attribute tags defined by the arguments {\tt WeightName1} and
1466	! {\tt WeightName2}, respectively. This paired integral is implicitly a
1467	! distributed operation (the whole motivation for pairing the integrals is
1468	! to reduce communication latency costs), and the Fortran MPI communicator
1469	! handle is defined by the input {\tt INTEGER} argument {\tt comm}. The
1470	! summation is an AllReduce operation, with all processes receiving the
1471	! global sum.
1472	!
1473	! {\bf N.B.: } The local lengths of the {\tt AttrVect} argument {\tt inAv1}
1474	! and the input {\tt REAL} array {\tt Weights1} must be equal. That is, there
1475	! must be a one-to-one correspondence between the field point values stored
1476	! in {\tt inAv1} and the point weights stored in {\tt Weights}. The same
1477	! relationship must apply between {\tt inAv2} and {\tt Weights2}.
1478	!
1479	! {\bf N.B.: } If {\tt SpatialIntegralRAttrV\_()} is invoked with the
1480	! optional {\tt LOGICAL} input argument {\tt SumWeights} set as {\tt .TRUE.},
1481	! then the value of {\tt WeightName1} must not conflict with any of the
1482	! {\tt REAL} attribute tags in {\tt inAv1} and the value of {\tt WeightName2}
1483	! must not conflict with any of the {\tt REAL} attribute tags in {\tt inAv2}.
1484	!
1485	! {\bf N.B.: } The output {\tt AttrVect} arguments {\tt outAv1} and
1486	! {\tt outAv2} are allocated data structures. The user must deallocate them
1487	! using the routine {\tt AttrVect\_clean()} when they are no longer needed.
1488	! Failure to do so will result in a memory leak.
1489	!
1490	! !INTERFACE:
1491
1492	subroutine PairedSpatialIntegralRAttrVSP_(inAv1, outAv1, Weights1, WeightName1, &
1493	inAv2, outAv2, Weights2, WeightName2, &
1494	SumWeights, comm)
1495	! ! USES:
1496
1497	use m_stdio
1498	use m_die
1499	use m_mpif90
1500	use m_realkinds, only : SP, FP
1501
1502	use m_AttrVect, only : AttrVect
1503	use m_AttrVect, only : AttrVect_lsize => lsize
1504	use m_AttrVect, only : AttrVect_nRAttr => nRAttr
1505
1506	use m_AttrVectReduce, only : AttrVect_LocalWeightedSumRAttr => &
1507	LocalWeightedSumRAttr
1508
1509	implicit none
1510
1511	! !INPUT PARAMETERS:
1512
1513	type(AttrVect), intent(IN) :: inAv1
1514	real(SP),dimension(:),pointer :: Weights1
1515	character(len=*), intent(IN) :: WeightName1
1516	type(AttrVect), intent(IN) :: inAv2
1517	real(SP),dimension(:),pointer :: Weights2
1518	character(len=*), intent(IN) :: WeightName2
1519	logical, optional, intent(IN) :: SumWeights
1520	integer, intent(IN) :: comm
1521
1522	! !OUTPUT PARAMETERS:
1523
1524	type(AttrVect), intent(OUT) :: outAv1
1525	type(AttrVect), intent(OUT) :: outAv2
1526
1527	! !REVISION HISTORY:
1528	! 10Jun02 - J.W. Larson <larson@mcs.anl.gov> - Initial version.
1529	!
1530	!EOP ___________________________________________________________________
1531
1532	character(len=*),parameter :: myname_=myname//'::PairedSpatialIntegralRAttrVSP_'
1533
1534	! Argument Sanity Checks:
1535
1536	integer :: ierr, length1, length2, PairedBufferLength
1537	integer :: nRA1, nRA2
1538	logical :: mySumWeights
1539	real(FP), dimension(:), pointer :: PairedBuffer, OutPairedBuffer
1540
1541	! Argument Validity Checks
1542
1543	if(AttrVect_lsize(inAv1) /= size(Weights1)) then
1544	ierr = AttrVect_lsize(inAv1) - size(Weights1)
1545	write(stderr,'(3a,i8,a,i8)') myname_, &
1546	':: inAv1 / Weights1 length mismatch: ', &
1547	' AttrVect_lsize(inAv1) = ',AttrVect_lsize(inAv1), &
1548	' size(Weights1) = ',size(Weights1)
1549	call die(myname_)
1550	endif
1551
1552	if(AttrVect_lsize(inAv2) /= size(Weights2)) then
1553	ierr = AttrVect_lsize(inAv2) - size(Weights2)
1554	write(stderr,'(3a,i8,a,i8)') myname_, &
1555	':: inAv2 / Weights2 length mismatch: ', &
1556	' AttrVect_lsize(inAv2) = ',AttrVect_lsize(inAv2), &
1557	' size(Weights2) = ',size(Weights2)
1558	call die(myname_)
1559	endif
1560
1561	! Are we summing the integration weights?
1562
1563	if(present(SumWeights)) then
1564	mySumWeights = SumWeights
1565	else
1566	mySumWeights = .FALSE.
1567	endif
1568
1569	! Compute the local contributions to the two integrals:
1570
1571	if(mySumWeights) then
1572	call AttrVect_LocalWeightedSumRAttr(inAv1, outAv1, Weights1, WeightName1)
1573	call AttrVect_LocalWeightedSumRAttr(inAv2, outAv2, Weights2, WeightName2)
1574	else
1575	call AttrVect_LocalWeightedSumRAttr(inAv1, outAv1, Weights1)
1576	call AttrVect_LocalWeightedSumRAttr(inAv2, outAv2, Weights2)
1577	endif
1578
1579	! Create the paired buffer for the Global Sum
1580
1581	nRA1 = AttrVect_nRAttr(outAv1)
1582	nRA2 = AttrVect_nRAttr(outAv2)
1583
1584	PairedBufferLength = nRA1 + nRA2
1585	allocate(PairedBuffer(PairedBufferLength), OutPairedBuffer(PairedBufferLength), &
1586	stat=ierr)
1587	if(ierr /= 0) then
1588	write(stderr,'(2a,i8)') myname_, &
1589	':: Fatal error--allocate(PairedBuffer...failed, ierr = ',ierr
1590	call die(myname_)
1591	endif
1592
1593	! Load the paired buffer
1594
1595	PairedBuffer(1:nRA1) = outAv1%rAttr(1:nRA1,1)
1596	PairedBuffer(nRA1+1:PairedBufferLength) = outAv2%rAttr(1:nRA2,1)
1597
1598	! Perform the global sum on the paired buffer
1599
1600	call MPI_AllReduce(PairedBuffer, OutPairedBuffer, PairedBufferLength, &
1601	MP_Type(PairedBuffer(1)), MP_SUM, comm, ierr)
1602	if(ierr /= 0) then
1603	write(stderr,'(2a,i8)') myname_, &
1604	':: Fatal Error--MPI_ALLREDUCE() failed with ierror = ',ierr
1605	call MP_perr_die(myname_,'MPI_ALLREDUCE() failed',ierr)
1606	endif
1607
1608	! Unload OutPairedBuffer into outAv1 and outAv2:
1609
1610	outAv1%rAttr(1:nRA1,1) = OutPairedBuffer(1:nRA1)
1611	outAv2%rAttr(1:nRA2,1) = OutPairedBuffer(nRA1+1:PairedBufferLength)
1612
1613	! Clean up allocated arrays:
1614
1615	deallocate(PairedBuffer, OutPairedBuffer, stat=ierr)
1616	if(ierr /= 0) then
1617	write(stderr,'(2a,i8)') myname_, &
1618	'ERROR--deallocate(PairedBuffer,...) failed, ierr = ',ierr
1619	call die(myname_)
1620	endif
1621
1622	end subroutine PairedSpatialIntegralRAttrVSP_
1623
1624	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1625	! Math and Computer Science Division, Argonne National Laboratory !
1626	! -------------------------------------------------------------------
1627	!
1628	! !IROUTINE: PairedSpatialIntegralRAttrVDP_ - Two spatial integrals.
1629	!
1630	! !DESCRIPTION:
1631	! Double precision interface version of PairedSpatialIntegralRAttrVSP_.
1632	!
1633	! !INTERFACE:
1634
1635	subroutine PairedSpatialIntegralRAttrVDP_(inAv1, outAv1, Weights1, WeightName1, &
1636	inAv2, outAv2, Weights2, WeightName2, &
1637	SumWeights, comm)
1638	! ! USES:
1639
1640	use m_stdio
1641	use m_die
1642	use m_mpif90
1643	use m_realkinds, only : DP, FP
1644
1645	use m_AttrVect, only : AttrVect
1646	use m_AttrVect, only : AttrVect_lsize => lsize
1647	use m_AttrVect, only : AttrVect_nRAttr => nRAttr
1648
1649	use m_AttrVectReduce, only : AttrVect_LocalWeightedSumRAttr => &
1650	LocalWeightedSumRAttr
1651
1652	implicit none
1653
1654	! !INPUT PARAMETERS:
1655
1656	type(AttrVect), intent(IN) :: inAv1
1657	real(DP),dimension(:),pointer :: Weights1
1658	character(len=*), intent(IN) :: WeightName1
1659	type(AttrVect), intent(IN) :: inAv2
1660	real(DP),dimension(:),pointer :: Weights2
1661	character(len=*), intent(IN) :: WeightName2
1662	logical, optional, intent(IN) :: SumWeights
1663	integer, intent(IN) :: comm
1664
1665	! !OUTPUT PARAMETERS:
1666
1667	type(AttrVect), intent(OUT) :: outAv1
1668	type(AttrVect), intent(OUT) :: outAv2
1669
1670	! !REVISION HISTORY:
1671	! 10Jun02 - J.W. Larson <larson@mcs.anl.gov> - Initial version.
1672	!
1673	! ___________________________________________________________________
1674
1675	character(len=*),parameter :: myname_=myname//'::PairedSpatialIntegralRAttrVDP_'
1676
1677	! Argument Sanity Checks:
1678
1679	integer :: ierr, length1, length2, PairedBufferLength
1680	integer :: nRA1, nRA2
1681	logical :: mySumWeights
1682	real(FP), dimension(:), pointer :: PairedBuffer, OutPairedBuffer
1683
1684	! Argument Validity Checks
1685
1686	if(AttrVect_lsize(inAv1) /= size(Weights1)) then
1687	ierr = AttrVect_lsize(inAv1) - size(Weights1)
1688	write(stderr,'(3a,i8,a,i8)') myname_, &
1689	':: inAv1 / Weights1 length mismatch: ', &
1690	' AttrVect_lsize(inAv1) = ',AttrVect_lsize(inAv1), &
1691	' size(Weights1) = ',size(Weights1)
1692	call die(myname_)
1693	endif
1694
1695	if(AttrVect_lsize(inAv2) /= size(Weights2)) then
1696	ierr = AttrVect_lsize(inAv2) - size(Weights2)
1697	write(stderr,'(3a,i8,a,i8)') myname_, &
1698	':: inAv2 / Weights2 length mismatch: ', &
1699	' AttrVect_lsize(inAv2) = ',AttrVect_lsize(inAv2), &
1700	' size(Weights2) = ',size(Weights2)
1701	call die(myname_)
1702	endif
1703
1704	! Are we summing the integration weights?
1705
1706	if(present(SumWeights)) then
1707	mySumWeights = SumWeights
1708	else
1709	mySumWeights = .FALSE.
1710	endif
1711
1712	! Compute the local contributions to the two integrals:
1713
1714	if(mySumWeights) then
1715	call AttrVect_LocalWeightedSumRAttr(inAv1, outAv1, Weights1, WeightName1)
1716	call AttrVect_LocalWeightedSumRAttr(inAv2, outAv2, Weights2, WeightName2)
1717	else
1718	call AttrVect_LocalWeightedSumRAttr(inAv1, outAv1, Weights1)
1719	call AttrVect_LocalWeightedSumRAttr(inAv2, outAv2, Weights2)
1720	endif
1721
1722	! Create the paired buffer for the Global Sum
1723
1724	nRA1 = AttrVect_nRAttr(outAv1)
1725	nRA2 = AttrVect_nRAttr(outAv2)
1726
1727	PairedBufferLength = nRA1 + nRA2
1728	allocate(PairedBuffer(PairedBufferLength), OutPairedBuffer(PairedBufferLength), &
1729	stat=ierr)
1730	if(ierr /= 0) then
1731	write(stderr,'(2a,i8)') myname_, &
1732	':: Fatal error--allocate(PairedBuffer...failed, ierr = ',ierr
1733	call die(myname_)
1734	endif
1735
1736	! Load the paired buffer
1737
1738	PairedBuffer(1:nRA1) = outAv1%rAttr(1:nRA1,1)
1739	PairedBuffer(nRA1+1:PairedBufferLength) = outAv2%rAttr(1:nRA2,1)
1740
1741	! Perform the global sum on the paired buffer
1742
1743	call MPI_AllReduce(PairedBuffer, OutPairedBuffer, PairedBufferLength, &
1744	MP_Type(PairedBuffer(1)), MP_SUM, comm, ierr)
1745	if(ierr /= 0) then
1746	write(stderr,'(2a,i8)') myname_, &
1747	':: Fatal Error--MPI_ALLREDUCE() failed with ierror = ',ierr
1748	call MP_perr_die(myname_,'MPI_ALLREDUCE() failed',ierr)
1749	endif
1750
1751	! Unload OutPairedBuffer into outAv1 and outAv2:
1752
1753	outAv1%rAttr(1:nRA1,1) = OutPairedBuffer(1:nRA1)
1754	outAv2%rAttr(1:nRA2,1) = OutPairedBuffer(nRA1+1:PairedBufferLength)
1755
1756	! Clean up allocated arrays:
1757
1758	deallocate(PairedBuffer, OutPairedBuffer, stat=ierr)
1759	if(ierr /= 0) then
1760	write(stderr,'(2a,i8)') myname_, &
1761	'ERROR--deallocate(PairedBuffer,...) failed, ierr = ',ierr
1762	call die(myname_)
1763	endif
1764
1765	end subroutine PairedSpatialIntegralRAttrVDP_
1766
1767	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1768	! Math and Computer Science Division, Argonne National Laboratory !
1769	!BOP -------------------------------------------------------------------
1770	!
1771	! !IROUTINE: PairedSpatialAverageRAttrVSP_ - Do two spatial averages at once.
1772	!
1773	! !DESCRIPTION:
1774	! This routine computes spatial averages of the {\tt REAL} attributes
1775	! of the {\tt REAL} attributes of the input {\tt AttrVect} arguments
1776	! {\tt inAv1} and {\tt inAv2}, returning the integrals in the output
1777	! {\tt AttrVect} arguments {\tt outAv1} and {\tt outAv2}, respectively .
1778	! The averages of {\tt inAv1} and {\tt inAv2} are computed using
1779	! spatial weights stored in the input {\tt REAL} array arguments
1780	! {\tt Weights1} and {\tt Weights2}, respectively. This paired average
1781	! is implicitly a
1782	! distributed operation (the whole motivation for pairing the integrals is
1783	! to reduce communication latency costs), and the Fortran MPI communicator
1784	! handle is defined by the input {\tt INTEGER} argument {\tt comm}. The
1785	! summation is an AllReduce operation, with all processes receiving the
1786	! global sum.
1787	!
1788	! {\bf N.B.: } The local lengths of the {\tt AttrVect} argument {\tt inAv1}
1789	! and the array {\tt Weights} must be equal. That is, there must be a
1790	! one-to-one correspondence between the field point values stored
1791	! in {\tt inAv1} and the spatial weights stored in {\tt Weights}
1792	!
1793	! {\bf N.B.: } The output {\tt AttrVect} arguments {\tt outAv1} and
1794	! {\tt outAv2} are allocated data structures. The user must deallocate them
1795	! using the routine {\tt AttrVect\_clean()} when they are no longer needed.
1796	! Failure to do so will result in a memory leak.
1797	!
1798	! !INTERFACE:
1799
1800	subroutine PairedSpatialAverageRAttrVSP_(inAv1, outAv1, Weights1, inAv2, &
1801	outAv2, Weights2, comm)
1802	! ! USES:
1803
1804	use m_stdio
1805	use m_die
1806	use m_mpif90
1807	use m_realkinds, only : SP, FP
1808
1809	use m_AttrVect, only : AttrVect
1810	use m_AttrVect, only : AttrVect_init => init
1811	use m_AttrVect, only : AttrVect_zero => zero
1812	use m_AttrVect, only : AttrVect_clean => clean
1813	use m_AttrVect, only : AttrVect_lsize => lsize
1814	use m_AttrVect, only : AttrVect_nRAttr => nRAttr
1815	use m_AttrVect, only : AttrVect_indexRA => indexRA
1816
1817	use m_AttrVectReduce, only : AttrVect_LocalWeightedSumRAttr => &
1818	LocalWeightedSumRAttr
1819
1820	use m_List, only : List
1821	use m_List, only : List_nullify => nullify
1822
1823	implicit none
1824
1825	! !INPUT PARAMETERS:
1826
1827	type(AttrVect), intent(IN) :: inAv1
1828	real(SP),dimension(:),pointer :: Weights1
1829	type(AttrVect), intent(IN) :: inAv2
1830	real(SP),dimension(:),pointer :: Weights2
1831	integer, intent(IN) :: comm
1832
1833	! !OUTPUT PARAMETERS:
1834
1835	type(AttrVect), intent(OUT) :: outAv1
1836	type(AttrVect), intent(OUT) :: outAv2
1837
1838	! !REVISION HISTORY:
1839	! 09May02 - J.W. Larson <larson@mcs.anl.gov> - Initial version.
1840	!
1841	!EOP ___________________________________________________________________
1842
1843	character(len=*),parameter :: myname_=myname//'::PairedSpatialAverageRAttrVSP_'
1844
1845	type(AttrVect) :: integratedAv1, integratedAv2
1846	type(List) :: nullIList
1847	integer :: i, ierr, iweight1, iweight2
1848
1849	! weight tags used to keep track of spatial weight sums
1850	character*8, parameter :: WeightName1='WeightSum1'
1851	character*8, parameter :: WeightName2='WeightSum2'
1852
1853	! Compute the paired spatial integral, including spatial weights:
1854
1855	call PairedSpatialIntegralsV(inAv1, integratedAv1, Weights1, WeightName1, &
1856	inAv2, integratedAv2, Weights2, WeightName2, &
1857	.TRUE., comm)
1858
1859	! Check value of summed weights (to avoid division by zero):
1860
1861	iweight1 = AttrVect_indexRA(integratedAv1, WeightName1)
1862	if(integratedAv1%rAttr(iweight1, 1) == 0._FP) then
1863	write(stderr,'(2a)') myname_, &
1864	'::ERROR--Global sum of grid weights in first integral is zero.'
1865	call die(myname_)
1866	endif
1867
1868	iweight2 = AttrVect_indexRA(integratedAv2, WeightName2)
1869	if(integratedAv2%rAttr(iweight2, 1) == 0._FP) then
1870	write(stderr,'(2a)') myname_, &
1871	'::ERROR--Global sum of grid weights in second integral is zero.'
1872	call die(myname_)
1873	endif
1874
1875	! Initialize output AttrVects outAv1 and outAv2:
1876
1877	call List_nullify(nullIList)
1878	call AttrVect_init(outAv1, iList=nullIList, rList=inAv1%rList, lsize=1)
1879	call AttrVect_zero(outAv1)
1880	call AttrVect_init(outAv2, iList=nullIList, rList=inAv2%rList, lsize=1)
1881	call AttrVect_zero(outAv2)
1882
1883	! Divide by global weight sum to compute spatial averages from
1884	! spatial integrals.
1885
1886	do i=1,AttrVect_nRAttr(outAv1)
1887	outAv1%rAttr(i,1) = integratedAv1%rAttr(i,1) &
1888	/ integratedAv1%rAttr(iweight1,1)
1889	end do
1890
1891	do i=1,AttrVect_nRAttr(outAv2)
1892	outAv2%rAttr(i,1) = integratedAv2%rAttr(i,1) &
1893	/ integratedAv2%rAttr(iweight2,1)
1894	end do
1895
1896	! Clean up temporary AttrVects:
1897
1898	call AttrVect_clean(integratedAv1)
1899	call AttrVect_clean(integratedAv2)
1900
1901	end subroutine PairedSpatialAverageRAttrVSP_
1902
1903	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1904	! Math and Computer Science Division, Argonne National Laboratory !
1905	! ----------------------------------------------------------------------
1906	!
1907	! !IROUTINE: PairedSpatialAverageRAttrVDP_ - Two spatial averages.
1908	!
1909	! !DESCRIPTION:
1910	! Double precision version of PairedSpatialAverageRAttrVSP_
1911	!
1912	! !INTERFACE:
1913
1914	subroutine PairedSpatialAverageRAttrVDP_(inAv1, outAv1, Weights1, inAv2, &
1915	outAv2, Weights2, comm)
1916	! ! USES:
1917
1918	use m_stdio
1919	use m_die
1920	use m_mpif90
1921	use m_realkinds, only : DP, FP
1922
1923	use m_AttrVect, only : AttrVect
1924	use m_AttrVect, only : AttrVect_init => init
1925	use m_AttrVect, only : AttrVect_zero => zero
1926	use m_AttrVect, only : AttrVect_clean => clean
1927	use m_AttrVect, only : AttrVect_lsize => lsize
1928	use m_AttrVect, only : AttrVect_nRAttr => nRAttr
1929	use m_AttrVect, only : AttrVect_indexRA => indexRA
1930
1931	use m_AttrVectReduce, only : AttrVect_LocalWeightedSumRAttr => &
1932	LocalWeightedSumRAttr
1933
1934	use m_List, only : List
1935	use m_List, only : List_nullify => nullify
1936
1937	implicit none
1938
1939	! !INPUT PARAMETERS:
1940
1941	type(AttrVect), intent(IN) :: inAv1
1942	real(DP),dimension(:),pointer :: Weights1
1943	type(AttrVect), intent(IN) :: inAv2
1944	real(DP),dimension(:),pointer :: Weights2
1945	integer, intent(IN) :: comm
1946
1947	! !OUTPUT PARAMETERS:
1948
1949	type(AttrVect), intent(OUT) :: outAv1
1950	type(AttrVect), intent(OUT) :: outAv2
1951
1952	! !REVISION HISTORY:
1953	! 09May02 - J.W. Larson <larson@mcs.anl.gov> - Initial version.
1954	!
1955	! ______________________________________________________________________
1956
1957	character(len=*),parameter :: myname_=myname//'::PairedSpatialAverageRAttrVDP_'
1958
1959	type(AttrVect) :: integratedAv1, integratedAv2
1960	type(List) :: nullIList
1961	integer :: i, ierr, iweight1, iweight2
1962
1963	! weight tags used to keep track of spatial weight sums
1964	character*8, parameter :: WeightName1='WeightSum1'
1965	character*8, parameter :: WeightName2='WeightSum2'
1966
1967	! Compute the paired spatial integral, including spatial weights:
1968
1969	call PairedSpatialIntegralsV(inAv1, integratedAv1, Weights1, WeightName1, &
1970	inAv2, integratedAv2, Weights2, WeightName2, &
1971	.TRUE., comm)
1972
1973	! Check value of summed weights (to avoid division by zero):
1974
1975	iweight1 = AttrVect_indexRA(integratedAv1, WeightName1)
1976	if(integratedAv1%rAttr(iweight1, 1) == 0._FP) then
1977	write(stderr,'(2a)') myname_, &
1978	'::ERROR--Global sum of grid weights in first integral is zero.'
1979	call die(myname_)
1980	endif
1981
1982	iweight2 = AttrVect_indexRA(integratedAv2, WeightName2)
1983	if(integratedAv2%rAttr(iweight2, 1) == 0._FP) then
1984	write(stderr,'(2a)') myname_, &
1985	'::ERROR--Global sum of grid weights in second integral is zero.'
1986	call die(myname_)
1987	endif
1988
1989	! Initialize output AttrVects outAv1 and outAv2:
1990
1991	call List_nullify(nullIList)
1992	call AttrVect_init(outAv1, iList=nullIList, rList=inAv1%rList, lsize=1)
1993	call AttrVect_zero(outAv1)
1994	call AttrVect_init(outAv2, iList=nullIList, rList=inAv2%rList, lsize=1)
1995	call AttrVect_zero(outAv2)
1996
1997	! Divide by global weight sum to compute spatial averages from
1998	! spatial integrals.
1999
2000	do i=1,AttrVect_nRAttr(outAv1)
2001	outAv1%rAttr(i,1) = integratedAv1%rAttr(i,1) &
2002	/ integratedAv1%rAttr(iweight1,1)
2003	end do
2004
2005	do i=1,AttrVect_nRAttr(outAv2)
2006	outAv2%rAttr(i,1) = integratedAv2%rAttr(i,1) &
2007	/ integratedAv2%rAttr(iweight2,1)
2008	end do
2009
2010	! Clean up temporary AttrVects:
2011
2012	call AttrVect_clean(integratedAv1)
2013	call AttrVect_clean(integratedAv2)
2014
2015	end subroutine PairedSpatialAverageRAttrVDP_
2016
2017	end module m_SpatialIntegralV

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: CPL/oasis3-mct/branches/OASIS3-MCT_2.0_branch/lib/mct/mct/m_SpatialIntegralV.F90 @ 4775

Download in other formats: