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