View Source

Gather/Scatter interfaces

Goal

Create generic interfaces for gather/scatter operations.

Motivation

Gathers and scatters between global fields and fields on the local physics or dynamics decomposition are done for

history output
restart I/O
reading initial files
reading boundary datasets

diagnostic output to log file
Lack of a generic interface leads to code like the following in the write_restart_history method of the history module:

select case (tape(t)%hlist(f)%field%decomp_type)
case (phys_decomp)
   call gather_chunk_to_field_hbuf (1, numlev, 1, plon, tape(t)%hlist(f)%hbuf, hbuf)
   call gather_chunk_to_field_int (1, 1, 1, plon, tape(t)%hlist(f)%nacs, fullnacs)
case (dyn_decomp)
   if ( dycore_is('LR') )then
# if ( defined STAGGERED )
h3.  NEW LR CODING
      lenarr = plon[numlev]plat
      if (tape(t)%hlist(f)%hbuf_prec == 8) then
         select case (numlev)
         case (1)
            call fv_gather('2d', tape(t)%hlist(f)%hbuf%buf8, lenarr, 2, hbuf%buf8)
         case (plev)
            call fv_gather('3dxzy', tape(t)%hlist(f)%hbuf%buf8, lenarr, 3, hbuf%buf8)
         case (plevp)
            call fv_gather('3dxzyp', tape(t)%hlist(f)%hbuf%buf8, lenarr, 3, hbuf%buf8)
         case default
            write(6,*)'WRITE_RESTART_HISTORY: invalid number of levels=', numlev
            call endrun ()
         end select
      else
         select case (numlev)
         case (1)
            call fv_gather4('2d', tape(t)%hlist(f)%hbuf%buf4, lenarr, 2, hbuf%buf4)
         case (plev)
            call fv_gather4('3dxzy', tape(t)%hlist(f)%hbuf%buf4, lenarr, 3, hbuf%buf4)
         case (plevp)
            call fv_gather4('3dxzyp', tape(t)%hlist(f)%hbuf%buf4, lenarr, 3, hbuf%buf4)
         case default
            write(6,*)'WRITE_RESTART_HISTORY: invalid number of levels=', numlev
            call endrun ()
         end select
      endif
      call fv_gatheri('2d', tape(t)%hlist(f)%nacs, lenarr, 2, fullnacs)
# endif
   else
      numowned = coldimin*numlev
      call compute_gsfactors (numowned, numsend, numrecv, displs)
      call mpigatherv_hbuf (tape(t)%hlist(f)%hbuf, numsend, mpireal, hbuf, numrecv, &
                            displs, mpireal, 0, mpicom)
      numowned = coldimin
      call compute_gsfactors (numowned, numsend, numrecv, displs)
      call mpigatherv (tape(t)%hlist(f)%nacs, numsend, mpiint, fullnacs, numrecv, &
                       displs, mpiint, 0, mpicom)
   endif
end select

All of this code is to gather two fields: hbuf and nacs. Similarly large sections of code are used to gather hbuf when the history file is written, and to scatter hbuf and nacs when the restart file is read. The dycore dependence of the gather/scatter methods means that code must be added to all three sections to support a new dycore.

Current interfaces for gathering physics chunks

Gathering physics chunks to global fields is done by the gather_chunk_to_hbuf methods in the phys_grid module:

   subroutine gather_chunk_to_field(fdim,mdim,ldim, &
                                     nlond,localchunks,globalfield)
   integer, intent(in) :: fdim      ! declared length of first dimension
   integer, intent(in) :: mdim      ! declared length of middle dimension
   integer, intent(in) :: ldim      ! declared length of last dimension
   integer, intent(in) :: nlond     ! declared number of longitudes
   real(r8), intent(in):: localchunks(fdim,pcols,mdim, &
                                      begchunk:endchunk,ldim) 
                                    ! local chunks
   real(r8), intent(out) :: globalfield(fdim,nlond,mdim,plat,ldim) 
                                    ! global field

This interface contains references to pcols, plat, begchunk, and endchunk which are present as module data in phys_grid. There are also real(r4) and integer versions which could easily be encapsulated in a generic interface.

Current interfaces for gathering FV dycore blocks

The fv_gather methods are contained in the io_dist module. They are implemented using the pargather method of the parutilities module.

   subroutine fv_gather(decomp_type, arr, lenarr, ndim, bufres)
      character(len=*) :: decomp_type
# if defined( SPMD )
      real(r8) arr(*)            ! Array to be gathered
# else
      real(r8) arr(lenarr)       ! Array (SMP-only)
# endif
      integer lenarr             ! Global length of array
      integer ndim               ! dimensionality (2 or 3) of array
      real(r8), intent(out) :: bufres(*)

Current interfaces for gathering spectral dycore blocks

The spectral dycores only implement a 1D decomposition by latitude slices. There is no single gather routine, but instead the gather is done by first calling the compute_gsfactors routine which sets up the indexing required by the mpigather methods. There is a wrapper routine mpigatherv_hbuf which encapsulates whether mpigatherv is called with the real(r8) or real(r4) version of hbuf. The compute_gsfactors method is contained in the spmd_dyn module and the mpigatherv_hbuf method is contained in the history module

   subroutine compute_gsfactors (numperlat, numtot, numperproc, displs)
      integer, intent(in)  :: numperlat            ! number of elements per latitude
      integer, intent(out) :: numtot               ! total number of elements (to send or recv)
      integer, intent(out) :: numperproc(0:npes-1) ! per-PE number of items to receive
      integer, intent(out) :: displs(0:npes-1)     ! per-PE displacements

   subroutine mpigatherv_hbuf (hbuf_send, numsend, mpireal1, hbuf_recv, &
                               numrecv, displs, mpireal2, root, comm)
      type (hbuffer_3d), intent(in   ) :: hbuf_send  ! send buffer
      type (hbuffer_3d), intent(inout) :: hbuf_recv  ! receive buffer
      integer :: numsend                ! number of items to be sent
      integer :: mpireal1               ! MPI real data type for hbuf_send
      integer :: mpireal2               ! MPI real data type for hbuf_recv
      integer :: numrecv(*)             ! number of items to be received
      integer :: displs(*)              ! displacement array
      integer, intent(in) :: root
      integer, intent(in) :: comm