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 17:57:11 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    18239232 :   integer function cam_filemap_getIndex(this)
     294             :     ! Dummy variable
     295             :     class(cam_filemap_t)                 :: this
     296             : 
     297    18239232 :     cam_filemap_getIndex = this%index
     298             : 
     299    18239232 :   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       38400 :   integer function cam_filemap_size(this)
     486             :     ! Dummy variable
     487             :     class(cam_filemap_t)                 :: this
     488             : 
     489       38400 :     if (associated(this%map)) then
     490       38400 :       cam_filemap_size = size(this%map,2)
     491             :     else
     492             :       cam_filemap_size = 0
     493             :     end if
     494       38400 :   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       24576 :   integer function cam_filemap_numMapped(this)
     524             :     ! Dummy variable
     525             :     class(cam_filemap_t)                 :: this
     526             : 
     527       24576 :     if (associated(this%map)) then
     528       24576 :       cam_filemap_numMapped = this%nmapped
     529             :     else
     530             :       cam_filemap_numMapped = 0
     531             :     end if
     532       24576 :   end function cam_filemap_numMapped
     533             : 
     534             :   !---------------------------------------------------------------------------
     535             :   !
     536             :   !  cam_filemap_mapVal: Calculate an offset value from a map
     537             :   !
     538             :   !---------------------------------------------------------------------------
     539    71958888 :   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    71958888 :     if (associated(this%map)) then
     551    71958888 :       if (present(dest_in)) then
     552    14140152 :         ndest = size(dest_in)
     553    14140152 :         d = 0
     554    42420456 :         d(:ndest) = dest_in(:ndest)
     555             :       else
     556   173456208 :         d = this%dest
     557             :       end if
     558    71958888 :       if (this%map(3, index) > 0) then
     559    67311196 :         cam_filemap_mapVal = ((this%map(3, index) - 1) * dsize(d(1))) + 1
     560    67311196 :         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    71958888 :   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       24576 :   subroutine cam_filemap_get_filemap(this, fieldlens, filelens, filemap,      &
     649       24576 :        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       24576 :     integer(iMap), allocatable    :: dsize(:)
     665             :     integer                       :: mapind(max_srcs)   ! Source index for map
     666       24576 :     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       24576 :     if (associated(filemap)) then
     679           0 :       deallocate(filemap)
     680             :       nullify(filemap)
     681             :     end if
     682             : 
     683             :     !
     684       66048 :     fileSize = product( int(filelens,kind=iMap) )
     685       24576 :     srccnt = size(fieldlens)
     686      112896 :     srclens(1:srccnt) = fieldlens(1:srccnt)
     687       24576 :     if (srccnt < 7) then
     688      157440 :       srclens(srccnt+1:7) = 1
     689             :     end if
     690             : 
     691             :     ! Allocate output filemap (dof)
     692      137472 :     allocate(filemap(product(fieldlens)))
     693    71983464 :     filemap = 0
     694             : 
     695             :     ! Find map source dimensions in input array
     696       24576 :     mapPos = 1                                ! To compare against map size
     697       24576 :     mapind = 0
     698       24576 :     if (present(src_in)) then
     699        3072 :       mapcnt = size(src_in)  ! Just used until end of loop below
     700        3072 :       if (mapcnt > max_srcs) then
     701           0 :         call endrun(subname//'src_in too large')
     702             :       end if
     703             :     end if
     704       73728 :     do j = 1, max_srcs
     705       49152 :       if (present(src_in)) then
     706        6144 :         if (mapcnt >= j) then
     707        6144 :           if (src_in(j) < 0) then
     708        3072 :             mapind(j) = srccnt + src_in(j) + 1
     709             :           else
     710        3072 :             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       43008 :         if (this%src(j) /= 0) then
     716             :           ! src < 0 is count from last array dim
     717       43008 :           if (this%src(j) < 0) then
     718       21504 :             mapind(j) = srccnt + this%src(j) + 1
     719             :           else
     720       21504 :             mapind(j) = this%src(j)
     721             :           end if
     722             :         end if
     723             :       end if
     724       73728 :       if (mapind(j) > 0) then
     725       49152 :         mapPos = mapPos * srclens(mapind(j))  ! To compare against map size
     726             :       end if
     727             :     end do
     728       73728 :     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       24576 :     mapSize = this%num_elem()
     734       24576 :     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       73728 :     allocate(dsize(size(filelens)))
     748       24576 :     dsize(1) = 1
     749       41472 :     do j = 1, size(filelens) - 1
     750       41472 :       dsize(j + 1) = dsize(j) * filelens(j)
     751             :     end do
     752             : 
     753             :     ! src_ind maps each source dimension to the corresponding file dimension
     754       73728 :     allocate(src_ind(srccnt))
     755       24576 :     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       88320 :       src_ind = 0
     762             :       fmind = 1      ! Here fmind is first available permutation slot
     763       88320 :       do j = 1, srccnt
     764             :         ! We need to find the offset location for each non-mapped src dimension
     765      142080 :         if (ANY(mapind == j)) then
     766             :           continue
     767             :         else
     768             :           ! Find the next available output dimension
     769       14592 :           if (present(dest_in)) then
     770        9216 :             do while (ANY(dest_in == fmind))
     771        3072 :               fmind = fmind + 1
     772        6144 :               if (fmind > size(dsize)) then
     773           0 :                 call endrun(subname//'permutation calculation dest_in error')
     774             :               end if
     775             :             end do
     776             :           else
     777       46080 :             do while (ANY(this%dest == fmind))
     778       11520 :               fmind = fmind + 1
     779       23040 :               if (fmind > size(dsize)) then
     780           0 :                 call endrun(subname//'permutation calculation dest error')
     781             :               end if
     782             :             end do
     783             :           end if
     784       14592 :           if (fmind > size(dsize)) then
     785           0 :             call endrun(subname//'permutation calculation error')
     786             :           end if
     787       14592 :           src_ind(j) = fmind
     788       14592 :           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       24576 :     fmind = 1
     795       49152 :     do i7 = 1, srclens(7)
     796       24576 :       i(7) = i7
     797       73728 :       do i6 = 1, srclens(6)
     798       24576 :         i(6) = i6
     799       73728 :         do i5 = 1, srclens(5)
     800       24576 :           i(5) = i5
     801       73728 :           do i4 = 1, srclens(4)
     802       24576 :             i(4) = i4
     803      127176 :             do i3 = 1, srclens(3)
     804       78024 :               i(3) = i3
     805     5259168 :               do i2 = 1, srclens(2)
     806     5156568 :                 i(2) = i2
     807    77193480 :                 do i1 = 1, srclens(1)
     808    71958888 :                   i(1) = i1
     809    71958888 :                   pos = 0  ! Offset from map pos so starts at zero
     810    71958888 :                   tind = 1 ! Offset into the map
     811    71958888 :                   tlen = 1
     812   287007264 :                   do j = 1, srccnt
     813   501227352 :                     if (ANY(mapind == j)) then
     814             :                       ! j is a distributed field dimension index
     815   143917776 :                       tind = tind + ((i(j) - 1) * tlen)
     816   143917776 :                       tlen = tlen * srclens(j)
     817             :                     else
     818             :                       ! j is a local field dimension index
     819    71130600 :                       pos = pos + ((i(j) - 1) * dsize(src_ind(j)))
     820             :                     end if
     821             :                   end do
     822    71958888 :                   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   129777624 :                   mapPos = this%map_val(tind, dsize, dest_in)
     829    71958888 :                   if ((mapPos > 0) .and. ((pos + mapPos) > fileSize)) then
     830           0 :                     call endrun(subname//'internal error, pos')
     831             :                   end if
     832    71958888 :                   if ((pos + mapPos) < 0) then
     833           0 :                     call endrun(subname//'internal error, mpos')
     834             :                   end if
     835    71958888 :                   if (mapPos > 0) then
     836    67311196 :                     filemap(fmind) = pos + mapPos
     837             :                   else
     838             :                     ! This element is not mapped
     839     4647692 :                     filemap(fmind) = 0
     840             :                   end if
     841    77115456 :                   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       24576 :     if ((fmind - 1) /= size(filemap)) then
     850           0 :       call endrun(subname//'internal error, fmind')
     851             :     end if
     852       24576 :     deallocate(dsize)
     853       24576 :   end subroutine cam_filemap_get_filemap
     854             : 
     855   585236880 :   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   585236880 :     if (.not. this%has_blocksize()) then
     862           0 :       call endrun('cam_filemap_get_blocksize: filemap has no blocks')
     863   585236880 :     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   585236880 :       cam_filemap_get_blocksize = this%blcksz(block_id)
     868             :     end if
     869   585236880 :   end function cam_filemap_get_blocksize
     870             : 
     871   585236880 :   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   585236880 :     cam_filemap_has_blocksize = associated(this%blcksz)
     879   585236880 :     if (present(lbnd)) then
     880           0 :       lbnd = LBOUND(this%blcksz, 1)
     881             :     end if
     882   585236880 :     if (present(ubnd)) then
     883           0 :       ubnd = UBOUND(this%blcksz, 1)
     884             :     end if
     885             : 
     886   585236880 :   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     2187264 :   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     2187264 :     rank = size(dims,1)
     907     2187264 :     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    15310848 :     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     6561792 :     do i = 1, max_srcs
     915     6561792 :       if (this%src(i) > 0) then
     916     2187264 :         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     2187264 :           if (associated(this%map)) then
     922     6561792 :             if (size(this%map) > 0) then
     923   147321552 :               dims(i, 1) = MINVAL(this%map(i,:))
     924   147321552 :               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     2187264 :       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     2187264 :         if (associated(this%map)) then
     938     6561792 :           if (size(this%map) > 0) then
     939   147321552 :             dims(rank + this%src(i) + 1, 1) = MINVAL(this%map(i,:))
     940   147321552 :             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     2187264 :   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