Line data Source code
1 : module cam_map_utils
2 : use pio, only: iMap=>PIO_OFFSET_KIND
3 : use cam_abortutils, only: endrun
4 :
5 : implicit none
6 : private
7 :
8 : public iMap
9 :
10 : integer, private, save :: unique_map_index = 0
11 : integer, private, parameter :: max_srcs = 2
12 : integer, private, parameter :: max_dests = 2
13 :
14 : !---------------------------------------------------------------------------
15 : !
16 : ! cam_filemap_t: Information for a CAM map (between data array and
17 : ! NetCDF file)
18 : !
19 : ! The map targets one or two dimensions of the NetCDF file.
20 : ! The 2-D map is useful for blocks of data (2 dimensions of array which
21 : ! map to one or two dimensions in the NetCDF file). For a 1-1 mapping,
22 : ! the first dimension is of size 3 (instead of 4).
23 : ! map(1, i) = index value of first src (e.g., lon, ncol)
24 : ! map(2, i) = index value of second src (e.g., lat, chunk)
25 : ! map(3, i) = global offset of first dest (e.g., lon, ncol)
26 : ! map(4, i) = global offset of second dest (e.g., lat, NA)
27 : ! src is the array dimension position corresponding to the map dimension
28 : ! in normal (not permuted) arrays.
29 : ! A negative pos denotes counting backwards from the last array dimension
30 : ! dest entries are the NetCDF dimension positions (must be increasing)
31 : ! It is an error to have src(1) == src(2).
32 : !
33 : !---------------------------------------------------------------------------
34 : type, public :: cam_filemap_t
35 : private
36 : integer :: index = 0
37 : integer(iMap), pointer :: map(:,:) => NULL()
38 : integer :: src(max_srcs) = 0
39 : integer :: dest(max_dests) = 0
40 : integer(iMap) :: limits(max_srcs,2) = -1
41 : integer :: nmapped = -1
42 : integer, pointer :: blcksz(:) => NULL() !e.g.,ncol(lchnk)
43 : contains
44 : procedure :: init => cam_filemap_init
45 : procedure :: get_index => cam_filemap_getIndex
46 : procedure :: new_index => cam_filemap_newIndex
47 : procedure :: clear => cam_filemap_clear
48 : procedure :: copy => cam_filemap_copy
49 : procedure :: copy_elem => cam_filemap_copyElem
50 : procedure :: num_elem => cam_filemap_size
51 : procedure :: is_mapped => cam_filemap_isMapped
52 : procedure :: num_mapped => cam_filemap_numMapped
53 : procedure :: map_val => cam_filemap_mapVal
54 : procedure :: coord_vals => cam_filemap_coordVals
55 : procedure :: coord_dests => cam_filemap_coordDests
56 : procedure :: get_filemap => cam_filemap_get_filemap
57 : procedure :: has_blocksize => cam_filemap_has_blocksize
58 : procedure :: blocksize => cam_filemap_get_blocksize
59 : procedure :: array_bounds => cam_filemap_get_array_bounds
60 : procedure :: active_cols => cam_filemap_get_active_cols
61 : procedure :: columnize => cam_filemap_columnize
62 : procedure :: compact => cam_filemap_compact
63 : !!XXgoldyXX: Cleanup when working
64 : ! procedure :: init_latlon => cam_filemap_init_latlon
65 : ! procedure :: init_unstruct => cam_filemap_init_unstruct
66 : ! generic, public :: init => init_latlon, init_unstruct
67 : end type cam_filemap_t
68 :
69 : !---------------------------------------------------------------------------
70 : !
71 : ! END: types BEGIN: private interfaces
72 : !
73 : !---------------------------------------------------------------------------
74 :
75 : contains
76 :
77 : !!#######################################################################
78 : !!
79 : !! index sorting routines:
80 : !! XXgoldyXX: Move to generic location?
81 : !!
82 : !!#######################################################################
83 :
84 0 : subroutine index_sort_vector(data, indices, compressval, dups_ok_in)
85 : use spmd_utils, only: mpi_integer, mpi_integer8, iam, mpicom, npes
86 : use shr_mpi_mod, only: shr_mpi_chkerr
87 : use cam_abortutils, only: endrun
88 : use m_MergeSorts, only: IndexSet, IndexSort
89 :
90 : ! Dummy arguments
91 : integer(iMap), pointer, intent(in) :: data(:)
92 : integer, pointer, intent(inout) :: indices(:)
93 : integer(iMap), optional, intent(in) :: compressval
94 : logical, optional, intent(in) :: dups_ok_in
95 :
96 : ! Local variables
97 : integer :: num_elem ! # mapped elements
98 : integer :: num_active ! # mapped pes
99 : integer :: ierr
100 : integer :: i
101 : integer :: lb, ub
102 : integer :: mycnt, my_first_elem
103 : integer :: mpi_group_world
104 : integer :: mpi_sort_group
105 : integer :: mpi_sort_comm
106 : integer :: my_sort_rank
107 0 : integer, allocatable :: sort_pes(:)
108 0 : integer, allocatable :: displs(:)
109 0 : integer, allocatable :: recvcounts(:)
110 0 : integer(iMap), allocatable :: elements(:)
111 0 : integer(iMap), allocatable :: local_elem(:)
112 0 : integer, allocatable :: ind(:)
113 0 : integer, allocatable :: temp(:)
114 : logical :: dups_ok ! .true. iff duplicates OK
115 : character(len=*), parameter :: subname = 'INDEX_SORT_VECTOR'
116 :
117 : ! Allow duplicate values?
118 0 : if (present(dups_ok_in)) then
119 0 : dups_ok = dups_ok_in
120 0 : if ((.not. dups_ok) .and. (.not. present(compressval))) then
121 0 : call endrun(trim(subname)//': dups_ok=.false. requires a compressval')
122 : end if
123 : else
124 : dups_ok = .true.
125 : end if
126 :
127 : ! The patch mapped values are in the number space of the master grid.
128 : ! They need to be compressed to go from 1 to the number of elements in the
129 : ! patch.
130 : ! Figure out the mapped elements in my portion of the patch mask
131 0 : if (.not. associated(data)) then
132 0 : mycnt = 0
133 0 : allocate(local_elem(mycnt))
134 0 : allocate(ind(mycnt))
135 0 : num_elem = 0
136 0 : lb = 0
137 0 : ub = -1
138 0 : else if (present(compressval)) then
139 0 : mycnt = COUNT(data /= compressval)
140 0 : allocate(local_elem(mycnt))
141 0 : allocate(ind(mycnt))
142 0 : num_elem = 0
143 0 : lb = LBOUND(data, 1)
144 0 : ub = UBOUND(data, 1)
145 0 : do i = lb, ub
146 0 : if (data(i) /= compressval) then
147 0 : num_elem = num_elem + 1
148 0 : local_elem(num_elem) = data(i)
149 0 : ind(num_elem) = i
150 : end if
151 : end do
152 : else
153 0 : lb = LBOUND(data, 1)
154 0 : ub = UBOUND(data, 1)
155 0 : mycnt = size(data)
156 0 : allocate(local_elem(mycnt))
157 0 : local_elem(1:mycnt) = data(lb:ub)
158 0 : num_elem = mycnt
159 : end if
160 :
161 : ! Find the tasks which have elements in this patch
162 : ! temp used for # elements per PE
163 0 : allocate(temp(0:npes-1))
164 : call MPI_allgather(mycnt, 1, MPI_integer, temp, 1, MPI_integer, &
165 0 : mpicom, ierr)
166 0 : call shr_mpi_chkerr(ierr, subname//': MPI_allgather elements')
167 0 : num_active = COUNT(temp > 0)
168 0 : if (num_active > 1) then
169 0 : allocate(sort_pes(num_active))
170 0 : allocate(recvcounts(num_active))
171 0 : allocate(displs(num_active))
172 0 : num_elem = 0
173 0 : displs(1) = 0
174 : ! Find the number of mapped elements and number of pes in this patch
175 0 : my_sort_rank = -1
176 0 : do i = 0, npes - 1
177 0 : if (temp(i) > 0) then
178 0 : num_elem = num_elem + 1
179 0 : if (num_elem > num_active) then
180 0 : call endrun(subname//": overrun of sort_pes array")
181 : end if
182 0 : sort_pes(num_elem) = i
183 0 : if (iam == i) then
184 0 : my_sort_rank = num_elem - 1
185 0 : my_first_elem = displs(num_elem) + 1
186 : end if
187 0 : recvcounts(num_elem) = temp(i)
188 0 : if (num_elem < num_active) then
189 0 : displs(num_elem + 1) = displs(num_elem) + recvcounts(num_elem)
190 : end if
191 : end if
192 : end do
193 0 : if (num_elem < num_active) then
194 0 : call endrun(subname//": underrun of sort_pes array")
195 : end if
196 0 : if (my_sort_rank >= 0) then
197 0 : num_elem = SUM(temp) ! Total number of elements to sort
198 : else
199 0 : num_elem = 0
200 : end if
201 0 : deallocate(temp) ! Cleanup
202 : ! Make a group with the active PEs
203 0 : call MPI_comm_group(mpicom, mpi_group_world, ierr)
204 0 : call shr_mpi_chkerr(ierr, subname//': MPI_comm_group mpi_group_world')
205 : call MPI_group_incl(mpi_group_world, num_active, sort_pes, &
206 0 : mpi_sort_group, ierr)
207 0 : call shr_mpi_chkerr(ierr, subname//': MPI_group_incl sort_pes')
208 : ! Make a new communicator with the active PEs
209 0 : call MPI_comm_create(mpicom, mpi_sort_group, mpi_sort_comm, ierr)
210 0 : call shr_mpi_chkerr(ierr, subname//': MPI_comm_create mpi_sort_comm')
211 : ! Collect all the elements for sorting (only active tasks now)
212 0 : allocate(elements(num_elem))
213 0 : if (mycnt > 0) then
214 : call MPI_allgatherv(local_elem, mycnt, MPI_integer8, &
215 0 : elements, recvcounts, displs, MPI_integer8, mpi_sort_comm, ierr)
216 0 : call shr_mpi_chkerr(ierr, subname//': MPI_allgatherv')
217 : ! Clean up for active PEs only
218 0 : call MPI_comm_free(mpi_sort_comm, ierr)
219 : end if
220 : ! General clean up
221 0 : call MPI_group_free(mpi_sort_group, ierr)
222 0 : deallocate(recvcounts)
223 0 : deallocate(displs)
224 0 : else if (mycnt > 0) then
225 : ! We are the only PE with patch info
226 0 : num_elem = mycnt
227 0 : allocate(elements(size(local_elem)))
228 0 : elements = local_elem
229 : my_first_elem = 1
230 : end if
231 : ! At this point, num_elem should always be the local # of items to sort
232 0 : if (num_elem > 0) then
233 : ! Sanity check
234 0 : if (size(elements) < num_elem) then
235 0 : call endrun(trim(subname)//": size(elements) must be >= num_elem")
236 : end if
237 : !! Do the sort, recvcounts will be the temporary index array
238 0 : allocate(recvcounts(size(elements)))
239 0 : allocate(displs(size(recvcounts)))
240 0 : call IndexSet(num_elem, recvcounts)
241 0 : call IndexSort(num_elem, recvcounts, elements, descend=.false.)
242 : ! Compress recvcounts (repeat data values)
243 0 : displs = 0
244 0 : do i = 1, num_elem - 1
245 0 : if (elements(recvcounts(i)) == elements(recvcounts(i + 1))) then
246 0 : displs(i + 1) = displs(i) + 1
247 : else
248 0 : displs(i + 1) = displs(i)
249 : end if
250 : end do
251 : ! Unload recvcounts into indices. Assume indices is initialized w/ default
252 0 : do i = 1, num_elem
253 0 : if ( (recvcounts(i) >= my_first_elem) .and. &
254 0 : (recvcounts(i) < (my_first_elem + mycnt))) then
255 0 : if (.not. dups_ok) then
256 : ! Use our indirect access to set the correct indices location
257 : ! ind array already has lb offset included
258 : ! Eliminate duplicate values
259 0 : if ((i > 1) .and. (displs(i) > displs(MAX((i - 1),1)))) then
260 0 : indices(ind(recvcounts(i) - my_first_elem + 1)) = compressval
261 : else
262 0 : indices(ind(recvcounts(i) - my_first_elem + 1)) = i - displs(i)
263 : end if
264 0 : else if (allocated(ind)) then
265 0 : indices(ind(recvcounts(i) - my_first_elem + 1)) = i - displs(i)
266 : else
267 : ! recvcounts points directly at a local location
268 : ! NB: repeat data values all get same index
269 0 : indices(recvcounts(i) - my_first_elem + lb) = i - displs(i)
270 : end if
271 : end if
272 : end do
273 0 : deallocate(recvcounts)
274 0 : deallocate(displs)
275 : end if
276 :
277 0 : if (allocated(ind)) then
278 0 : deallocate(ind)
279 : end if
280 0 : if (allocated(elements)) then
281 0 : deallocate(elements)
282 : end if
283 0 : deallocate(local_elem)
284 :
285 0 : end subroutine index_sort_vector
286 :
287 : !!#######################################################################
288 : !!
289 : !! CAM grid mapping functions
290 : !!
291 : !!#######################################################################
292 :
293 10791936 : integer function cam_filemap_getIndex(this)
294 : ! Dummy variable
295 : class(cam_filemap_t) :: this
296 :
297 10791936 : cam_filemap_getIndex = this%index
298 :
299 10791936 : end function cam_filemap_getIndex
300 :
301 7680 : subroutine cam_filemap_newIndex(this)
302 : ! Dummy variable
303 : class(cam_filemap_t) :: this
304 :
305 7680 : unique_map_index = unique_map_index + 1
306 7680 : this%index = unique_map_index
307 :
308 7680 : end subroutine cam_filemap_newIndex
309 :
310 7680 : subroutine cam_filemap_init(this, pemap, unstruct, src, dest)
311 : ! Dummy arguments
312 : class(cam_filemap_t) :: this
313 : integer(iMap), pointer :: pemap(:,:) ! Map elem for this PE
314 : logical, intent(in) :: unstruct
315 : integer, intent(in) :: src(:)
316 : integer, optional, intent(in) :: dest(:)
317 :
318 : ! Local variables
319 : integer :: i ! Loop index
320 : integer :: index
321 :
322 : ! This shouldn't happen but maybe we will decide to reuse these
323 7680 : if (associated(this%map)) then
324 0 : deallocate(this%map)
325 0 : nullify(this%map)
326 : end if
327 :
328 : ! Check in case these ever change (because algorithm then must be modified)
329 : if ((max_srcs /= 2) .or. (max_dests /= 2)) then
330 : call endrun('cam_filemap_init: max_src or max_dest modified')
331 : end if
332 :
333 : ! Some items are simply copied
334 7680 : if (associated(pemap)) then
335 7680 : this%map => pemap
336 : else
337 0 : nullify(this%map)
338 : end if
339 30720 : this%src = src
340 7680 : if (present(dest)) then
341 : ! Structred grids will likely always have dest = (1, 2) but maybe . . .
342 30720 : this%dest = dest
343 0 : else if (unstruct) then
344 : ! Unstructured grids are (so far) always spread along the first dimension
345 0 : this%dest(1) = 1
346 0 : this%dest(2) = 0
347 : else
348 0 : this%dest(1) = 1
349 0 : this%dest(2) = 2
350 : end if
351 : ! We may have holes in the 'block' decomposition which is specified by
352 : ! having src(2) < 0.
353 : ! NB: This is currently a special purpose hack in that it is purely
354 : ! convention that the last dimension specifies the block index and
355 : ! that those blocks may not be filled.
356 : ! The proper way to generalize this functionality is to allow
357 : ! src(1) to also be < 0 and to look for holes there as well
358 7680 : if (associated(this%map)) then
359 23040 : do i = 1, max_srcs
360 23040 : if (ANY(this%map(i,:) > 0)) then
361 : ! Min of all src(i) values
362 1308864 : this%limits(i, 1) = MINVAL(this%map(i,:), mask=(this%map(i,:)>0))
363 : ! Max of all src(i) values
364 1293504 : this%limits(i, 2) = MAXVAL(this%map(i,:), mask=(this%map(i,:)>0))
365 : else
366 0 : this%limits(i,1) = 0
367 0 : this%limits(i,2) = -1
368 : end if
369 : end do
370 :
371 7680 : this%nmapped = 0
372 646752 : do index = 1, this%num_elem()
373 646752 : if ((this%dest(1) > 0) .and. (this%map(max_srcs+1, index) > 0)) then
374 486008 : if (this%dest(2) > 0) then
375 : ! Can't do this test unless we know the dim 2 is large enough
376 0 : if (this%map(max_srcs+2, index) > 0) then
377 0 : this%nmapped = this%nmapped + 1
378 : end if
379 : else
380 486008 : this%nmapped = this%nmapped + 1
381 : end if
382 : end if
383 : end do
384 : else
385 0 : this%limits(:,1) = 0
386 0 : this%limits(:,2) = -1
387 0 : this%nmapped = 0
388 : end if
389 7680 : if (src(max_srcs) < 0) then
390 : ! This shouldn't happen but maybe we will decide to reuse these
391 7680 : if (associated(this%blcksz)) then
392 0 : deallocate(this%blcksz)
393 0 : nullify(this%blcksz)
394 : end if
395 23040 : allocate(this%blcksz(this%limits(max_srcs,1):this%limits(max_srcs,2)))
396 57072 : this%blcksz = 0
397 646752 : do i = 1, this%num_elem()
398 639072 : index = this%map(max_srcs, i)
399 646752 : if (this%is_mapped(i)) then
400 486008 : this%blcksz(index) = this%blcksz(index) + 1
401 : end if
402 : end do
403 : end if
404 :
405 7680 : call this%new_index()
406 :
407 7680 : end subroutine cam_filemap_init
408 :
409 0 : subroutine cam_filemap_clear(this)
410 : ! Dummy arguments
411 : class(cam_filemap_t) :: this
412 :
413 0 : if (associated(this%map)) then
414 0 : this%map = 0
415 : end if
416 0 : this%limits(:,1) = 0
417 0 : this%limits(:,2) = -1
418 0 : this%nmapped = 0
419 :
420 : ! Update the index in case a decomp was made with the old values
421 0 : call this%new_index()
422 :
423 0 : end subroutine cam_filemap_clear
424 :
425 0 : subroutine cam_filemap_copy(this, map)
426 : ! Dummy arguments
427 : class(cam_filemap_t) :: this
428 : type(cam_filemap_t), intent(in) :: map
429 :
430 : ! This shouldn't happen but maybe we will decide to reuse these
431 0 : if (associated(this%map)) then
432 0 : deallocate(this%map)
433 0 : nullify(this%map)
434 : end if
435 :
436 0 : if (associated(map%map)) then
437 0 : allocate(this%map(size(map%map, 1), size(map%map, 2)))
438 0 : if (map%num_elem() > 0) then
439 0 : this%map = map%map
440 : end if
441 : else
442 0 : nullify(this%map)
443 : end if
444 0 : this%src = map%src
445 0 : this%dest = map%dest
446 0 : this%limits = map%limits
447 0 : this%nmapped = map%nmapped
448 :
449 : ! This shouldn't happen but maybe we will decide to reuse these
450 0 : if (associated(this%blcksz)) then
451 0 : deallocate(this%blcksz)
452 0 : nullify(this%blcksz)
453 : end if
454 0 : if (associated(map%blcksz)) then
455 0 : allocate(this%blcksz(size(map%blcksz)))
456 0 : this%blcksz = map%blcksz
457 : end if
458 :
459 : ! Even a copy has to have a unique index
460 0 : call this%new_index()
461 :
462 0 : end subroutine cam_filemap_copy
463 :
464 0 : subroutine cam_filemap_copyElem(this, map, index)
465 : ! Dummy arguments
466 : class(cam_filemap_t) :: this
467 : type(cam_filemap_t), intent(in) :: map
468 : integer, intent(in) :: index
469 :
470 0 : if (this%is_mapped(index)) then
471 0 : this%nmapped = this%nmapped - 1
472 : end if
473 0 : this%map(:, index) = map%map(:, index)
474 0 : if (this%is_mapped(index)) then
475 0 : this%nmapped = this%nmapped + 1
476 : end if
477 :
478 0 : end subroutine cam_filemap_copyElem
479 :
480 : !---------------------------------------------------------------------------
481 : !
482 : ! cam_filemap_size: Total number of elements in the map
483 : !
484 : !---------------------------------------------------------------------------
485 33792 : integer function cam_filemap_size(this)
486 : ! Dummy variable
487 : class(cam_filemap_t) :: this
488 :
489 33792 : if (associated(this%map)) then
490 33792 : cam_filemap_size = size(this%map,2)
491 : else
492 : cam_filemap_size = 0
493 : end if
494 33792 : end function cam_filemap_size
495 :
496 : !---------------------------------------------------------------------------
497 : !
498 : ! cam_filemap_isMapped: Return .true. iff the index value is mapped
499 : !
500 : !---------------------------------------------------------------------------
501 639072 : elemental logical function cam_filemap_isMapped(this, index)
502 : ! Dummy arguments
503 : class(cam_filemap_t), intent(in) :: this
504 : integer, intent(in) :: index
505 :
506 639072 : if (associated(this%map)) then
507 639072 : cam_filemap_isMapped = (this%map(3, index) > 0)
508 639072 : if ((size(this%map, 1) > 3) .and. cam_filemap_isMapped) then
509 0 : cam_filemap_isMapped = (this%map(4, index) > 0)
510 : ! No else needed
511 : end if
512 : else
513 : cam_filemap_isMapped = .false.
514 : end if
515 :
516 639072 : end function cam_filemap_isMapped
517 :
518 : !---------------------------------------------------------------------------
519 : !
520 : ! cam_filemap_numMapped: Number of elements in the map with non-zero entries
521 : !
522 : !---------------------------------------------------------------------------
523 18432 : integer function cam_filemap_numMapped(this)
524 : ! Dummy variable
525 : class(cam_filemap_t) :: this
526 :
527 18432 : if (associated(this%map)) then
528 18432 : cam_filemap_numMapped = this%nmapped
529 : else
530 : cam_filemap_numMapped = 0
531 : end if
532 18432 : end function cam_filemap_numMapped
533 :
534 : !---------------------------------------------------------------------------
535 : !
536 : ! cam_filemap_mapVal: Calculate an offset value from a map
537 : !
538 : !---------------------------------------------------------------------------
539 17528400 : integer(iMap) function cam_filemap_mapVal(this, index, dsize, dest_in)
540 : ! Dummy arguments
541 : class(cam_filemap_t) :: this
542 : integer, intent(in) :: index
543 : integer(iMap), intent(in) :: dsize(:)
544 : integer, optional, intent(in) :: dest_in(:)
545 :
546 : ! Local variables
547 : integer :: ndest
548 : integer :: d(max_dests)
549 :
550 17528400 : if (associated(this%map)) then
551 17528400 : if (present(dest_in)) then
552 3510000 : ndest = size(dest_in)
553 3510000 : d = 0
554 10530000 : d(:ndest) = dest_in(:ndest)
555 : else
556 42055200 : d = this%dest
557 : end if
558 17528400 : if (this%map(3, index) > 0) then
559 16183860 : cam_filemap_mapVal = ((this%map(3, index) - 1) * dsize(d(1))) + 1
560 16183860 : if (size(this%map, 1) > 3) then
561 0 : if (this%map(4, index) > 0) then
562 0 : cam_filemap_mapVal = ((this%map(4, index) - 1) * dsize(d(2))) + &
563 0 : cam_filemap_mapVal ! No +1 because it is offset from map(3,:)
564 : else
565 : cam_filemap_mapVal = 0
566 : end if
567 : ! No else needed
568 : end if
569 : else
570 : cam_filemap_mapVal = 0
571 : end if
572 : else
573 : cam_filemap_mapVal = 0
574 : end if
575 17528400 : end function cam_filemap_mapVal
576 :
577 : !---------------------------------------------------------------------------
578 : !
579 : ! cam_filemap_coordVals: Find coord indices matching map index
580 : !
581 : !---------------------------------------------------------------------------
582 0 : subroutine cam_filemap_coordVals(this, index, lonIndex, latIndex, isMapped)
583 : ! Dummy arguments
584 : class(cam_filemap_t) :: this
585 : integer, intent(in) :: index
586 : integer, intent(out) :: lonIndex
587 : integer, intent(out) :: latIndex
588 : logical, optional, intent(out) :: isMapped
589 :
590 0 : if (associated(this%map)) then
591 0 : if (size(this%map,1) > (max_srcs + 1)) then
592 0 : lonIndex = this%map(1, index)
593 0 : latIndex = this%map(2, index)
594 : else
595 0 : lonIndex = index
596 0 : latIndex = index
597 : end if
598 0 : if (present(isMapped)) then
599 0 : isMapped = this%is_mapped(index)
600 : end if
601 0 : else if (present(isMapped)) then
602 0 : lonIndex = 0
603 0 : latIndex = 0
604 0 : isMapped = .false.
605 : else
606 0 : call endrun("cam_filemap_coordVals: must have map or pass isMapped")
607 : end if
608 0 : end subroutine cam_filemap_coordVals
609 :
610 : !---------------------------------------------------------------------------
611 : !
612 : ! cam_filemap_coordDests: Find coord indices matching map index
613 : !
614 : !---------------------------------------------------------------------------
615 0 : subroutine cam_filemap_coordDests(this, index, lonIndex, latIndex, isMapped)
616 : ! Dummy arguments
617 : class(cam_filemap_t) :: this
618 : integer, intent(in) :: index
619 : integer, intent(out) :: lonIndex
620 : integer, intent(out) :: latIndex
621 : logical, optional, intent(out) :: isMapped
622 :
623 0 : if (associated(this%map)) then
624 0 : lonIndex = this%map(max_srcs + 1, index)
625 0 : if (size(this%map,1) > (max_srcs + 1)) then
626 0 : latIndex = this%map(max_srcs + 2, index)
627 : else
628 0 : latIndex = 0
629 : end if
630 0 : if (present(isMapped)) then
631 0 : isMapped = this%is_mapped(index)
632 : end if
633 0 : else if (present(isMapped)) then
634 0 : lonIndex = 0
635 0 : latIndex = 0
636 0 : isMapped = .false.
637 : else
638 0 : call endrun("cam_filemap_coordDests: must have map or pass isMapped")
639 : end if
640 0 : end subroutine cam_filemap_coordDests
641 :
642 : !---------------------------------------------------------------------------
643 : !
644 : ! cam_filemap_get_filemap: Create the mapping between an array and a file
645 : ! This is the workhorse function for creating a PIO decomp DOF
646 : !
647 : !---------------------------------------------------------------------------
648 18432 : subroutine cam_filemap_get_filemap(this, fieldlens, filelens, filemap, &
649 18432 : src_in, dest_in, permutation_in)
650 : use shr_kind_mod, only: SHR_KIND_CL
651 :
652 : ! Dummy arguments
653 : class(cam_filemap_t) :: this
654 : integer, intent(in) :: fieldlens(:)
655 : integer, intent(in) :: filelens(:)
656 : integer(iMap), pointer :: filemap(:)
657 : integer, optional, intent(in) :: src_in(:)
658 : integer, optional, intent(in) :: dest_in(:)
659 : integer, optional, intent(in) :: permutation_in(:)
660 :
661 : ! Local variables
662 : integer :: srclens(7) ! field dim lens
663 : integer :: srccnt ! Rank of fieldlens
664 18432 : integer(iMap), allocatable :: dsize(:)
665 : integer :: mapind(max_srcs) ! Source index for map
666 18432 : integer, allocatable :: src_ind(:) ! Map src to file dims
667 : integer :: fmind, j
668 : integer(iMap) :: mapSize, mapPos, pos, fileSize
669 : integer :: mapcnt ! Dimension count
670 : integer :: tind, tlen ! Temporarys
671 : integer :: i1, i2, i3, i4, i5, i6, i7
672 : integer :: i(7)
673 : character(len=SHR_KIND_CL) :: errmsg
674 : character(len=32) :: errfmt
675 : character(len=*), parameter :: subname = 'cam_filemap_get_filemap: '
676 :
677 : ! This shouldn't happen but, who knows what evil lurks in the hearts of SEs
678 18432 : if (associated(filemap)) then
679 0 : deallocate(filemap)
680 : nullify(filemap)
681 : end if
682 :
683 : !
684 50688 : fileSize = product( int(filelens,kind=iMap) )
685 18432 : srccnt = size(fieldlens)
686 85248 : srclens(1:srccnt) = fieldlens(1:srccnt)
687 18432 : if (srccnt < 7) then
688 117504 : srclens(srccnt+1:7) = 1
689 : end if
690 :
691 : ! Allocate output filemap (dof)
692 103680 : allocate(filemap(product(fieldlens)))
693 17546832 : filemap = 0
694 :
695 : ! Find map source dimensions in input array
696 18432 : mapPos = 1 ! To compare against map size
697 18432 : mapind = 0
698 18432 : if (present(src_in)) then
699 1536 : mapcnt = size(src_in) ! Just used until end of loop below
700 1536 : if (mapcnt > max_srcs) then
701 0 : call endrun(subname//'src_in too large')
702 : end if
703 : end if
704 55296 : do j = 1, max_srcs
705 36864 : if (present(src_in)) then
706 3072 : if (mapcnt >= j) then
707 3072 : if (src_in(j) < 0) then
708 1536 : mapind(j) = srccnt + src_in(j) + 1
709 : else
710 1536 : mapind(j) = src_in(j)
711 : end if
712 : end if
713 : else
714 : ! A src == 0 means that map dimension is not used
715 33792 : if (this%src(j) /= 0) then
716 : ! src < 0 is count from last array dim
717 33792 : if (this%src(j) < 0) then
718 16896 : mapind(j) = srccnt + this%src(j) + 1
719 : else
720 16896 : mapind(j) = this%src(j)
721 : end if
722 : end if
723 : end if
724 55296 : if (mapind(j) > 0) then
725 36864 : mapPos = mapPos * srclens(mapind(j)) ! To compare against map size
726 : end if
727 : end do
728 55296 : mapcnt = COUNT(mapind /= 0)
729 :
730 : ! Check that map matches dims
731 : ! Since patch maps are compressed, we can't do an equal compare but
732 : ! it is still an error if the map has more elements than the array
733 18432 : mapSize = this%num_elem()
734 18432 : if (mapPos < this%num_mapped()) then
735 0 : if (mapcnt > 1) then
736 0 : write(errfmt, '(a,i0,2a)') "(a,i0,a,", mapcnt, '(i0,", "),")")'
737 : else
738 0 : write(errfmt, '(a,i0,2a)') '(a,i0,a,i0,")")'
739 : end if
740 0 : write(errmsg, errfmt) 'Map size (', &
741 0 : this%num_mapped(), ') too large for array dims (', &
742 0 : srclens(mapind(1:mapcnt))
743 0 : call endrun(subname//trim(errmsg))
744 : end if
745 :
746 : ! dsize is a global offset for each dimension
747 55296 : allocate(dsize(size(filelens)))
748 18432 : dsize(1) = 1
749 32256 : do j = 1, size(filelens) - 1
750 32256 : dsize(j + 1) = dsize(j) * filelens(j)
751 : end do
752 :
753 : ! src_ind maps each source dimension to the corresponding file dimension
754 55296 : allocate(src_ind(srccnt))
755 18432 : if (present(permutation_in)) then
756 0 : if (size(permutation_in) /= size(src_ind)) then
757 0 : call endrun(subname//'permutation_in must have same rank as fieldlens')
758 : end if
759 0 : src_ind = permutation_in
760 : else
761 66816 : src_ind = 0
762 : fmind = 1 ! Here fmind is first available permutation slot
763 66816 : do j = 1, srccnt
764 : ! We need to find the offset location for each non-mapped src dimension
765 108288 : if (ANY(mapind == j)) then
766 : continue
767 : else
768 : ! Find the next available output dimension
769 11520 : if (present(dest_in)) then
770 4608 : do while (ANY(dest_in == fmind))
771 1536 : fmind = fmind + 1
772 3072 : if (fmind > size(dsize)) then
773 0 : call endrun(subname//'permutation calculation dest_in error')
774 : end if
775 : end do
776 : else
777 39936 : do while (ANY(this%dest == fmind))
778 9984 : fmind = fmind + 1
779 19968 : if (fmind > size(dsize)) then
780 0 : call endrun(subname//'permutation calculation dest error')
781 : end if
782 : end do
783 : end if
784 11520 : if (fmind > size(dsize)) then
785 0 : call endrun(subname//'permutation calculation error')
786 : end if
787 11520 : src_ind(j) = fmind
788 11520 : fmind = fmind + 1
789 : end if
790 : end do
791 : end if
792 :
793 : ! Step through the map and fill in local positions for each entry
794 18432 : fmind = 1
795 36864 : do i7 = 1, srclens(7)
796 18432 : i(7) = i7
797 55296 : do i6 = 1, srclens(6)
798 18432 : i(6) = i6
799 55296 : do i5 = 1, srclens(5)
800 18432 : i(5) = i5
801 55296 : do i4 = 1, srclens(4)
802 18432 : i(4) = i4
803 99432 : do i3 = 1, srclens(3)
804 62568 : i(3) = i3
805 1360800 : do i2 = 1, srclens(2)
806 1279800 : i(2) = i2
807 18870768 : do i1 = 1, srclens(1)
808 17528400 : i(1) = i1
809 17528400 : pos = 0 ! Offset from map pos so starts at zero
810 17528400 : tind = 1 ! Offset into the map
811 17528400 : tlen = 1
812 69520320 : do j = 1, srccnt
813 120918960 : if (ANY(mapind == j)) then
814 : ! j is a distributed field dimension index
815 35056800 : tind = tind + ((i(j) - 1) * tlen)
816 35056800 : tlen = tlen * srclens(j)
817 : else
818 : ! j is a local field dimension index
819 16935120 : pos = pos + ((i(j) - 1) * dsize(src_ind(j)))
820 : end if
821 : end do
822 17528400 : if (tind > mapSize) then
823 : write(errmsg, '(2(a,i0),a,12x,5(i0,", "),")")') &
824 0 : 'internal error, tind (', tind, ') > mapSize (', &
825 0 : mapSize, '), srclens = (', srclens(1:5)
826 0 : call endrun(subname//trim(errmsg))
827 : end if
828 31546800 : mapPos = this%map_val(tind, dsize, dest_in)
829 17528400 : if ((mapPos > 0) .and. ((pos + mapPos) > fileSize)) then
830 0 : call endrun(subname//'internal error, pos')
831 : end if
832 17528400 : if ((pos + mapPos) < 0) then
833 0 : call endrun(subname//'internal error, mpos')
834 : end if
835 17528400 : if (mapPos > 0) then
836 16183860 : filemap(fmind) = pos + mapPos
837 : else
838 : ! This element is not mapped
839 1344540 : filemap(fmind) = 0
840 : end if
841 18808200 : fmind = fmind + 1
842 : end do
843 : end do
844 : end do
845 : end do
846 : end do
847 : end do
848 : end do
849 18432 : if ((fmind - 1) /= size(filemap)) then
850 0 : call endrun(subname//'internal error, fmind')
851 : end if
852 18432 : deallocate(dsize)
853 18432 : end subroutine cam_filemap_get_filemap
854 :
855 231686064 : integer function cam_filemap_get_blocksize(this, block_id)
856 :
857 : ! Dummy arguments
858 : class(cam_filemap_t) :: this
859 : integer, intent(in) :: block_id
860 :
861 231686064 : if (.not. this%has_blocksize()) then
862 0 : call endrun('cam_filemap_get_blocksize: filemap has no blocks')
863 231686064 : else if ((block_id < LBOUND(this%blcksz, 1)) .or. &
864 : (block_id > UBOUND(this%blcksz, 1))) then
865 0 : call endrun('cam_filemap_get_blocksize: block_id out of range')
866 : else
867 231686064 : cam_filemap_get_blocksize = this%blcksz(block_id)
868 : end if
869 231686064 : end function cam_filemap_get_blocksize
870 :
871 231686064 : logical function cam_filemap_has_blocksize(this, lbnd, ubnd)
872 :
873 : ! Dummy arguments
874 : class(cam_filemap_t) :: this
875 : integer, optional, intent(out) :: lbnd
876 : integer, optional, intent(out) :: ubnd
877 :
878 231686064 : cam_filemap_has_blocksize = associated(this%blcksz)
879 231686064 : if (present(lbnd)) then
880 0 : lbnd = LBOUND(this%blcksz, 1)
881 : end if
882 231686064 : if (present(ubnd)) then
883 0 : ubnd = UBOUND(this%blcksz, 1)
884 : end if
885 :
886 231686064 : end function cam_filemap_has_blocksize
887 :
888 : !---------------------------------------------------------------------------
889 : !
890 : ! cam_filemap_get_array_bounds: Sets grid bounds for the relevant array
891 : ! Only modifies the dimensions corresponding to the map's src
892 : ! dims should be sized (rank,2) with the second dimension used
893 : ! to store lower(1) and upper(2) bounds
894 : !
895 : !---------------------------------------------------------------------------
896 1052928 : subroutine cam_filemap_get_array_bounds(this, dims)
897 :
898 : ! Dummy arguments
899 : class(cam_filemap_t) :: this
900 : integer, intent(inout) :: dims(:,:)
901 :
902 : ! Local variables
903 : integer :: rank ! rank of target array
904 : integer :: i ! Loop variable
905 :
906 1052928 : rank = size(dims,1)
907 1052928 : if (size(dims,2) < 2) then
908 0 : call endrun('cam_filemap_get_array_bounds: second dim of dims must be 2')
909 : end if
910 :
911 7370496 : if (MAXVAL(this%limits(1:max_srcs, 1:2)) > HUGE(kind(dims))) then
912 0 : call endrun('cam_filemap_get_array_bounds: limits too large')
913 : end if
914 3158784 : do i = 1, max_srcs
915 3158784 : if (this%src(i) > 0) then
916 1052928 : if (this%src(i) > rank) then
917 0 : call endrun('cam_filemap_get_array_bounds: rank too small')
918 : else
919 : !!XXgoldyXX: Maybe modify definition of this%limits?
920 : ! dims(i, 1:2) = INT(this%limits(i, 1:2), kind=kind(dims))
921 1052928 : if (associated(this%map)) then
922 3158784 : if (size(this%map) > 0) then
923 71519184 : dims(i, 1) = MINVAL(this%map(i,:))
924 71519184 : dims(i, 2) = MAXVAL(this%map(i,:))
925 : else
926 0 : dims(i, 1) = 0
927 0 : dims(i, 2) = -1
928 : end if
929 : else
930 0 : dims(i, 1) = 0
931 0 : dims(i, 2) = -1
932 : end if
933 : end if
934 1052928 : else if (this%src(i) < 0) then
935 : !!XXgoldyXX: Maybe modify definition of this%limits?
936 : ! dims(rank + this%src(i) + 1, 1:2) = INT(this%limits(i, 1:2), kind=kind(dims))
937 1052928 : if (associated(this%map)) then
938 3158784 : if (size(this%map) > 0) then
939 71519184 : dims(rank + this%src(i) + 1, 1) = MINVAL(this%map(i,:))
940 71519184 : dims(rank + this%src(i) + 1, 2) = MAXVAL(this%map(i,:))
941 : else
942 0 : dims(rank + this%src(i) + 1, 1) = 0
943 0 : dims(rank + this%src(i) + 1, 2) = -1
944 : end if
945 : else
946 0 : dims(rank + this%src(i) + 1, 1) = 0
947 0 : dims(rank + this%src(i) + 1, 2) = -1
948 : end if
949 : else
950 : ! src(i)==0 means unused position
951 0 : dims(i,:) = 0
952 : end if
953 : end do
954 1052928 : end subroutine cam_filemap_get_array_bounds
955 :
956 : !---------------------------------------------------------------------------
957 : !
958 : ! cam_filemap_get_active_cols: Find which columns are active in a dimension
959 : ! Because we normally decompose columns (blocks or chunks) in the
960 : ! last dimension, we default to that dimension
961 : !
962 : !---------------------------------------------------------------------------
963 0 : subroutine cam_filemap_get_active_cols(this, colnum, active, srcdim_in)
964 :
965 : ! Dummy arguments
966 : class(cam_filemap_t) :: this
967 : integer, intent(in) :: colnum
968 : logical, intent(out) :: active(:)
969 : integer, optional, intent(in) :: srcdim_in
970 :
971 : ! Local variables
972 : integer :: srcdim
973 :
974 0 : if (present(srcdim_in)) then
975 0 : srcdim = srcdim_in
976 : else
977 : srcdim = max_srcs
978 : end if
979 :
980 : ! Sanity checks
981 0 : if ((srcdim < 1) .or. (srcdim > max_srcs)) then
982 0 : call endrun('cam_filemap_get_active_cols: srcdim out of range')
983 0 : else if (this%src(srcdim) >= 0) then
984 0 : call endrun('cam_filemap_get_active_cols: Invalid srcdim')
985 0 : else if (size(active) < size(this%map, (3 - srcdim))) then
986 0 : call endrun('cam_filemap_get_active_cols: active too small')
987 0 : else if (colnum < LBOUND(this%map, srcdim)) then
988 0 : call endrun('cam_filemap_get_active_cols: colnum too small')
989 0 : else if (colnum > UBOUND(this%map, srcdim)) then
990 0 : call endrun('cam_filemap_get_active_cols: colnum too large')
991 : ! else OK
992 : end if
993 :
994 0 : active = .false.
995 : !!XXgoldyXX: This is probably completely wrong. What we want is column info
996 0 : select case(srcdim)
997 : case (1)
998 0 : active(1:size(this%map, 2)) = (this%map(colnum,:) > 0)
999 : case (2)
1000 0 : active(1:size(this%map, 1)) = (this%map(:,colnum) > 0)
1001 : case default
1002 0 : call endrun('cam_filemap_get_active_cols: Invalid srcdim?!?')
1003 : end select
1004 0 : end subroutine cam_filemap_get_active_cols
1005 :
1006 : !---------------------------------------------------------------------------
1007 : !
1008 : ! cam_filemap_columnize: Convert lon/lat map to ncol
1009 : !
1010 : !---------------------------------------------------------------------------
1011 0 : subroutine cam_filemap_columnize(this)
1012 : use spmd_utils, only: mpi_sum, mpi_integer, mpicom
1013 : use shr_mpi_mod, only: shr_mpi_chkerr
1014 :
1015 : ! Dummy argument
1016 : class(cam_filemap_t) :: this
1017 :
1018 : ! Local variables
1019 : integer :: i, j
1020 : integer :: lmax(2)
1021 : integer :: maxind(2)
1022 : integer :: offset(2)
1023 : integer :: ierr
1024 : integer(iMap), pointer :: newmap(:,:) => NULL()
1025 : character(len=*), parameter :: subname = 'CAM_FILEMAP_COLUMNIZE'
1026 : ! Create a new map with same size and ordering
1027 : ! We need the max lon/lat for calculating global offsets
1028 0 : lmax = 0
1029 0 : if (associated(this%map)) then
1030 0 : if (size(this%map, 1) == max_srcs) then
1031 0 : call endrun(trim(subname)//': must have at least 1 destination coord')
1032 0 : else if (size(this%map, 1) > (max_srcs + 2)) then
1033 0 : call endrun(trim(subname)//': has more than 2 destination coords')
1034 : end if
1035 0 : if (size(this%map, 2) > 0) then
1036 0 : lmax(1) = MAXVAL(this%map(max_srcs + 1, :))
1037 0 : if (size(this%map, 1) == (max_srcs + 2)) then
1038 0 : lmax(2) = MAXVAL(this%map(max_srcs + 2, :))
1039 : end if
1040 : end if
1041 : end if
1042 0 : call MPI_allreduce(lmax(1), maxind(1), 1, MPI_integer, mpi_sum, mpicom, ierr)
1043 0 : call shr_mpi_chkerr(ierr, subname//': MPI_allreduce maxlon')
1044 0 : call MPI_allreduce(lmax(2), maxind(2), 1, MPI_integer, mpi_sum, mpicom, ierr)
1045 0 : call shr_mpi_chkerr(ierr, subname//': MPI_allreduce maxlat')
1046 0 : if (associated(this%map)) then
1047 0 : if (size(this%map, 1) == (max_srcs + 2))then
1048 : ! Create the new map
1049 0 : allocate(newmap(max_srcs + 1, size(this%map, 2)))
1050 : ! Who's on first?
1051 0 : if (ANY(this%dest(1:2) <= 0)) then
1052 0 : call endrun(trim(subname)//': can only handle positive dest indices')
1053 0 : else if (this%dest(1) < this%dest(2)) then
1054 0 : offset(1) = 1
1055 0 : offset(2) = maxind(1)
1056 0 : else if (this%dest(1) > this%dest(2)) then
1057 0 : offset(1) = maxind(2)
1058 0 : offset(2) = 1
1059 : else
1060 0 : call endrun(trim(subname)//': dest indices cannot be equal')
1061 : end if
1062 0 : do i = 1, size(newmap, 2)
1063 0 : newmap(1:max_srcs, i) = this%map(1:max_srcs, i)
1064 0 : j = ((this%map(max_srcs+1, i) * offset(1)) + &
1065 0 : (this%map(max_srcs+2, i) * offset(2)))
1066 0 : newmap(max_srcs+1, i) = j
1067 : end do
1068 : ! Replace our map with the new one
1069 0 : deallocate(this%map)
1070 0 : this%map => newmap
1071 0 : nullify(newmap)
1072 : ! Fixup dest
1073 0 : this%dest(1) = 1
1074 0 : this%dest(2:) = 0
1075 : ! no else (do nothing if already ncol)
1076 : end if ! End if lon/lat map
1077 : end if
1078 :
1079 0 : end subroutine cam_filemap_columnize
1080 :
1081 : !---------------------------------------------------------------------------
1082 : !
1083 : ! cam_filemap_compact: Pack all active elements into a 1-based file order
1084 : ! Also pack latitude and longitude maps
1085 : !
1086 : !---------------------------------------------------------------------------
1087 0 : subroutine cam_filemap_compact(this, lonmap, latmap, &
1088 : num_lons, num_lats, num_mapped, columnize, dups_ok_in)
1089 : use spmd_utils, only: mpi_sum, mpi_integer, mpicom
1090 : use shr_mpi_mod, only: shr_mpi_chkerr
1091 :
1092 : ! Dummy arguments
1093 : class(cam_filemap_t) :: this
1094 : integer(iMap), pointer :: lonmap(:)
1095 : integer(iMap), pointer :: latmap(:)
1096 : integer, optional, intent(out) :: num_lons
1097 : integer, optional, intent(out) :: num_lats
1098 : integer, optional, intent(out) :: num_mapped
1099 : logical, optional, intent(in) :: columnize ! Convert to ncol
1100 : logical, optional, intent(in) :: dups_ok_in ! Dup coords OK
1101 :
1102 : ! Local variables
1103 : integer :: i
1104 : integer :: ierr
1105 : integer(iMap), pointer :: data(:) => NULL()
1106 : integer, pointer :: indices(:) => NULL()
1107 : integer :: tmp_size
1108 : logical :: dok
1109 : character(len=*), parameter :: subname = 'CAM_FILEMAP_COMPACT'
1110 :
1111 : !! Possibly convert lon/lat map to ncol
1112 0 : if (present(columnize)) then
1113 0 : if (columnize) then
1114 0 : call this%columnize()
1115 : end if
1116 : end if
1117 : !! Are duplicate coordinate indices (lat/lon) OK?
1118 0 : if (present(dups_ok_in)) then
1119 0 : dok = dups_ok_in
1120 : else
1121 0 : dok = .false.
1122 : end if
1123 : ! Get a global index sort of mapped elements.
1124 0 : do i = 1, max_dests
1125 0 : if (this%dest(i) > 0) then
1126 0 : if (associated(this%map)) then
1127 0 : if (size(this%map, 1) >= max_srcs + i) then
1128 0 : data => this%map(max_srcs + i, :)
1129 : else
1130 0 : nullify(data)
1131 : end if
1132 : else
1133 0 : nullify(data)
1134 : end if
1135 : ! Allocate indices if necessary
1136 0 : if (associated(indices)) then
1137 0 : deallocate(indices)
1138 : nullify(indices)
1139 : end if
1140 0 : if (associated(data)) then
1141 0 : allocate(indices(LBOUND(data, 1):UBOUND(data, 1)))
1142 0 : indices = 0
1143 : else
1144 0 : allocate(indices(0))
1145 : end if
1146 : end if
1147 0 : call index_sort_vector(data, indices, compressval=0_iMap)
1148 0 : if (associated(data) .and. associated(indices)) then
1149 0 : data = indices
1150 : end if
1151 : end do
1152 0 : if (associated(indices)) then
1153 0 : deallocate(indices)
1154 : nullify(indices)
1155 : end if
1156 : ! Get a global index sort of lat and lon maps
1157 : !! Compress latmap
1158 0 : if (associated(latmap)) then
1159 : ! Allocate indices
1160 0 : allocate(indices(LBOUND(latmap, 1):UBOUND(latmap, 1)))
1161 0 : indices = 0
1162 : end if
1163 0 : call index_sort_vector(latmap, indices, compressval=0_iMap, dups_ok_in=dok)
1164 0 : if (associated(latmap)) then
1165 0 : latmap = indices
1166 0 : deallocate(indices)
1167 : nullify(indices)
1168 : end if
1169 : !! Compress lonmap
1170 0 : if (associated(lonmap)) then
1171 0 : allocate(indices(LBOUND(lonmap, 1):UBOUND(lonmap, 1)))
1172 0 : indices = 0
1173 : end if
1174 0 : call index_sort_vector(lonmap, indices, compressval=0_iMap, dups_ok_in=dok)
1175 0 : if (associated(lonmap)) then
1176 0 : lonmap = indices
1177 0 : deallocate(indices)
1178 : nullify(indices)
1179 : end if
1180 :
1181 0 : if (present(num_mapped)) then
1182 0 : if (this%num_elem() > 0) then
1183 : ! Total number of mapped elements
1184 0 : allocate(indices(this%num_elem()))
1185 0 : indices = [ (i, i = 1, size(indices)) ]
1186 0 : tmp_size = COUNT(this%is_mapped(indices))
1187 0 : deallocate(indices)
1188 : nullify(indices)
1189 : else
1190 0 : tmp_size = 0
1191 : end if
1192 : call MPI_allreduce(tmp_size, num_mapped, 1, MPI_integer, &
1193 0 : mpi_sum, mpicom, ierr)
1194 0 : call shr_mpi_chkerr(ierr, subname//': MPI_allreduce num_mapped')
1195 0 : if (num_mapped <= 0) then
1196 0 : call endrun(trim(subname)//': num_mapped <= 0')
1197 : end if
1198 : end if
1199 0 : if (present(num_lons)) then
1200 0 : if (associated(lonmap)) then
1201 0 : tmp_size = COUNT(lonmap /= 0)
1202 : else
1203 0 : tmp_size = 0
1204 : end if
1205 : call MPI_allreduce(tmp_size, num_lons, 1, MPI_integer, &
1206 0 : mpi_sum, mpicom, ierr)
1207 0 : call shr_mpi_chkerr(ierr, subname//': MPI_allreduce num_lons')
1208 0 : if (num_lons <= 0) then
1209 0 : call endrun(trim(subname)//': numlons <= 0')
1210 : end if
1211 : end if
1212 0 : if (present(num_lats)) then
1213 0 : if (associated(latmap)) then
1214 0 : tmp_size = COUNT(latmap /= 0)
1215 : else
1216 0 : tmp_size = 0
1217 : end if
1218 : call MPI_allreduce(tmp_size, num_lats, 1, MPI_integer, &
1219 0 : mpi_sum, mpicom, ierr)
1220 0 : call shr_mpi_chkerr(ierr, subname//': MPI_allreduce num_lats')
1221 0 : if (num_lats <= 0) then
1222 0 : call endrun(trim(subname)//': numlats <= 0')
1223 : end if
1224 : end if
1225 :
1226 0 : end subroutine cam_filemap_compact
1227 :
1228 0 : end module cam_map_utils
|
