LCOV - code coverage report
Current view: | top level - utils - cam_map_utils.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 170 | 532 | 32.0 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 11 | 21 | 52.4 %

              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    43342848 :   integer function cam_filemap_getIndex(this)
         294             :     ! Dummy variable
         295             :     class(cam_filemap_t)                 :: this
         296             :
         297    43342848 :     cam_filemap_getIndex = this%index
         298             :
         299    43342848 :   end function cam_filemap_getIndex
         300             :
         301        6912 :   subroutine cam_filemap_newIndex(this)
         302             :     ! Dummy variable
         303             :     class(cam_filemap_t)                 :: this
         304             :
         305        6912 :     unique_map_index = unique_map_index + 1
         306        6912 :     this%index = unique_map_index
         307             :
         308        6912 :   end subroutine cam_filemap_newIndex
         309             :
         310        6912 :   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        6912 :     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        6912 :     if (associated(pemap)) then
         335        6912 :       this%map => pemap
         336             :     else
         337           0 :       nullify(this%map)
         338             :     end if
         339       27648 :     this%src =  src
         340        6912 :     if (present(dest)) then
         341             :       ! Structred grids will likely always have dest = (1, 2) but maybe . . .
         342       27648 :       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        6912 :     if (associated(this%map)) then
         359       20736 :       do i = 1, max_srcs
         360       20736 :         if (ANY(this%map(i,:) > 0)) then
         361             :           ! Min of all src(i) values
         362     1132992 :           this%limits(i, 1) = MINVAL(this%map(i,:), mask=(this%map(i,:)>0))
         363             :           ! Max of all src(i) values
         364     1119168 :           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        6912 :       this%nmapped = 0
         372      559584 :       do index = 1, this%num_elem()
         373      559584 :         if ((this%dest(1) > 0) .and. (this%map(max_srcs+1, index) > 0)) then
         374      437406 :           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      437406 :             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        6912 :     if (src(max_srcs) < 0) then
         390             :       ! This shouldn't happen but maybe we will decide to reuse these
         391        6912 :       if (associated(this%blcksz)) then
         392           0 :         deallocate(this%blcksz)
         393           0 :         nullify(this%blcksz)
         394             :       end if
         395       20736 :       allocate(this%blcksz(this%limits(max_srcs,1):this%limits(max_srcs,2)))
         396       50904 :       this%blcksz = 0
         397      559584 :       do i = 1, this%num_elem()
         398      552672 :         index = this%map(max_srcs, i)
         399      559584 :         if (this%is_mapped(i)) then
         400      437406 :           this%blcksz(index) = this%blcksz(index) + 1
         401             :         end if
         402             :       end do
         403             :     end if
         404             :
         405        6912 :     call this%new_index()
         406             :
         407        6912 :   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       46848 :   integer function cam_filemap_size(this)
         486             :     ! Dummy variable
         487             :     class(cam_filemap_t)                 :: this
         488             :
         489       46848 :     if (associated(this%map)) then
         490       46848 :       cam_filemap_size = size(this%map,2)
         491             :     else
         492             :       cam_filemap_size = 0
         493             :     end if
         494       46848 :   end function cam_filemap_size
         495             :
         496             :   !---------------------------------------------------------------------------
         497             :   !
         498             :   !  cam_filemap_isMapped: Return .true. iff the index value is mapped
         499             :   !
         500             :   !---------------------------------------------------------------------------
         501      552672 :   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      552672 :     if (associated(this%map)) then
         507      552672 :       cam_filemap_isMapped = (this%map(3, index) > 0)
         508      552672 :       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      552672 :   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       33024 :   integer function cam_filemap_numMapped(this)
         524             :     ! Dummy variable
         525             :     class(cam_filemap_t)                 :: this
         526             :
         527       33024 :     if (associated(this%map)) then
         528       33024 :       cam_filemap_numMapped = this%nmapped
         529             :     else
         530             :       cam_filemap_numMapped = 0
         531             :     end if
         532       33024 :   end function cam_filemap_numMapped
         533             :
         534             :   !---------------------------------------------------------------------------
         535             :   !
         536             :   !  cam_filemap_mapVal: Calculate an offset value from a map
         537             :   !
         538             :   !---------------------------------------------------------------------------
         539   138039912 :   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   138039912 :     if (associated(this%map)) then
         551   138039912 :       if (present(dest_in)) then
         552    20480760 :         ndest = size(dest_in)
         553    20480760 :         d = 0
         554    61442280 :         d(:ndest) = dest_in(:ndest)
         555             :       else
         556   352677456 :         d = this%dest
         557             :       end if
         558   138039912 :       if (this%map(3, index) > 0) then
         559   132143596 :         cam_filemap_mapVal = ((this%map(3, index) - 1) * dsize(d(1))) + 1
         560   132143596 :         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   138039912 :   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       33024 :   subroutine cam_filemap_get_filemap(this, fieldlens, filelens, filemap,      &
         649       33024 :        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       33024 :     integer(iMap), allocatable    :: dsize(:)
         665             :     integer                       :: mapind(max_srcs)   ! Source index for map
         666       33024 :     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       33024 :     if (associated(filemap)) then
         679           0 :       deallocate(filemap)
         680             :       nullify(filemap)
         681             :     end if
         682             :
         683             :     !
         684       91392 :     fileSize = product( int(filelens,kind=iMap) )
         685       33024 :     srccnt = size(fieldlens)
         686      155136 :     srclens(1:srccnt) = fieldlens(1:srccnt)
         687       33024 :     if (srccnt < 7) then
         688      208128 :       srclens(srccnt+1:7) = 1
         689             :     end if
         690             :
         691             :     ! Allocate output filemap (dof)
         692      188160 :     allocate(filemap(product(fieldlens)))
         693   138072936 :     filemap = 0
         694             :
         695             :     ! Find map source dimensions in input array
         696       33024 :     mapPos = 1                                ! To compare against map size
         697       33024 :     mapind = 0
         698       33024 :     if (present(src_in)) then
         699        9216 :       mapcnt = size(src_in)  ! Just used until end of loop below
         700        9216 :       if (mapcnt > max_srcs) then
         701           0 :         call endrun(subname//'src_in too large')
         702             :       end if
         703             :     end if
         704       99072 :     do j = 1, max_srcs
         705       66048 :       if (present(src_in)) then
         706       18432 :         if (mapcnt >= j) then
         707       18432 :           if (src_in(j) < 0) then
         708        9216 :             mapind(j) = srccnt + src_in(j) + 1
         709             :           else
         710        9216 :             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       47616 :         if (this%src(j) /= 0) then
         716             :           ! src < 0 is count from last array dim
         717       47616 :           if (this%src(j) < 0) then
         718       23808 :             mapind(j) = srccnt + this%src(j) + 1
         719             :           else
         720       23808 :             mapind(j) = this%src(j)
         721             :           end if
         722             :         end if
         723             :       end if
         724       99072 :       if (mapind(j) > 0) then
         725       66048 :         mapPos = mapPos * srclens(mapind(j))  ! To compare against map size
         726             :       end if
         727             :     end do
         728       99072 :     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       33024 :     mapSize = this%num_elem()
         734       33024 :     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       99072 :     allocate(dsize(size(filelens)))
         748       33024 :     dsize(1) = 1
         749       58368 :     do j = 1, size(filelens) - 1
         750       58368 :       dsize(j + 1) = dsize(j) * filelens(j)
         751             :     end do
         752             :
         753             :     ! src_ind maps each source dimension to the corresponding file dimension
         754       99072 :     allocate(src_ind(srccnt))
         755       33024 :     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      122112 :       src_ind = 0
         762             :       fmind = 1      ! Here fmind is first available permutation slot
         763      122112 :       do j = 1, srccnt
         764             :         ! We need to find the offset location for each non-mapped src dimension
         765      201216 :         if (ANY(mapind == j)) then
         766             :           continue
         767             :         else
         768             :           ! Find the next available output dimension
         769       23040 :           if (present(dest_in)) then
         770       27648 :             do while (ANY(dest_in == fmind))
         771        9216 :               fmind = fmind + 1
         772       18432 :               if (fmind > size(dsize)) then
         773           0 :                 call endrun(subname//'permutation calculation dest_in error')
         774             :               end if
         775             :             end do
         776             :           else
         777       55296 :             do while (ANY(this%dest == fmind))
         778       13824 :               fmind = fmind + 1
         779       27648 :               if (fmind > size(dsize)) then
         780           0 :                 call endrun(subname//'permutation calculation dest error')
         781             :               end if
         782             :             end do
         783             :           end if
         784       23040 :           if (fmind > size(dsize)) then
         785           0 :             call endrun(subname//'permutation calculation error')
         786             :           end if
         787       23040 :           src_ind(j) = fmind
         788       23040 :           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       33024 :     fmind = 1
         795       66048 :     do i7 = 1, srclens(7)
         796       33024 :       i(7) = i7
         797       99072 :       do i6 = 1, srclens(6)
         798       33024 :         i(6) = i6
         799       99072 :         do i5 = 1, srclens(5)
         800       33024 :           i(5) = i5
         801       99072 :           do i4 = 1, srclens(4)
         802       33024 :             i(4) = i4
         803      178128 :             do i3 = 1, srclens(3)
         804      112080 :               i(3) = i3
         805     9431736 :               do i2 = 1, srclens(2)
         806     9286632 :                 i(2) = i2
         807   147438624 :                 do i1 = 1, srclens(1)
         808   138039912 :                   i(1) = i1
         809   138039912 :                   pos = 0  ! Offset from map pos so starts at zero
         810   138039912 :                   tind = 1 ! Offset into the map
         811   138039912 :                   tlen = 1
         812   551331360 :                   do j = 1, srccnt
         813   963794520 :                     if (ANY(mapind == j)) then
         814             :                       ! j is a distributed field dimension index
         815   276079824 :                       tind = tind + ((i(j) - 1) * tlen)
         816   276079824 :                       tlen = tlen * srclens(j)
         817             :                     else
         818             :                       ! j is a local field dimension index
         819   137211624 :                       pos = pos + ((i(j) - 1) * dsize(src_ind(j)))
         820             :                     end if
         821             :                   end do
         822   138039912 :                   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   255599064 :                   mapPos = this%map_val(tind, dsize, dest_in)
         829   138039912 :                   if ((mapPos > 0) .and. ((pos + mapPos) > fileSize)) then
         830           0 :                     call endrun(subname//'internal error, pos')
         831             :                   end if
         832   138039912 :                   if ((pos + mapPos) < 0) then
         833           0 :                     call endrun(subname//'internal error, mpos')
         834             :                   end if
         835   138039912 :                   if (mapPos > 0) then
         836   132143596 :                     filemap(fmind) = pos + mapPos
         837             :                   else
         838             :                     ! This element is not mapped
         839     5896316 :                     filemap(fmind) = 0
         840             :                   end if
         841   147326544 :                   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       33024 :     if ((fmind - 1) /= size(filemap)) then
         850           0 :       call endrun(subname//'internal error, fmind')
         851             :     end if
         852       33024 :     deallocate(dsize)
         853       33024 :   end subroutine cam_filemap_get_filemap
         854             :
         855  1184517216 :   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  1184517216 :     if (.not. this%has_blocksize()) then
         862           0 :       call endrun('cam_filemap_get_blocksize: filemap has no blocks')
         863  1184517216 :     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  1184517216 :       cam_filemap_get_blocksize = this%blcksz(block_id)
         868             :     end if
         869  1184517216 :   end function cam_filemap_get_blocksize
         870             :
         871  1184517216 :   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  1184517216 :     cam_filemap_has_blocksize = associated(this%blcksz)
         879  1184517216 :     if (present(lbnd)) then
         880           0 :       lbnd = LBOUND(this%blcksz, 1)
         881             :     end if
         882  1184517216 :     if (present(ubnd)) then
         883           0 :       ubnd = UBOUND(this%blcksz, 1)
         884             :     end if
         885             :
         886  1184517216 :   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     5179392 :   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     5179392 :     rank = size(dims,1)
         907     5179392 :     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    36255744 :     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    15538176 :     do i = 1, max_srcs
         915    15538176 :       if (this%src(i) > 0) then
         916     5179392 :         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     5179392 :           if (associated(this%map)) then
         922    15538176 :             if (size(this%map) > 0) then
         923   346455216 :               dims(i, 1) = MINVAL(this%map(i,:))
         924   346455216 :               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     5179392 :       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     5179392 :         if (associated(this%map)) then
         938    15538176 :           if (size(this%map) > 0) then
         939   346455216 :             dims(rank + this%src(i) + 1, 1) = MINVAL(this%map(i,:))
         940   346455216 :             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     5179392 :   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

Generated by: LCOV version 1.14